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Abstract. Weakly interacting massive particles (WIMPs) are amongst the most interesting 
dark matter (DM) candidates. Many DM candidates naturally arise in theories beyond the 
standard model (SM) of particle physics, like weak-scale supersymmetry (SUSY). Experi- 
ments aim to detect WIMPs by scattering, annihilation or direct production, and thereby 
determine the underlying theory to which they belong, along with its parameters. Here we 
examine the prospects for further constraining the Constrained Minimal Supersymmetric 
Standard Model (CMSSM) with future ton-scale direct detection experiments. We consider 
ton-scale extrapolations of three current experiments: CDMS, XENON and COUPP, with 
1000 kg-years of raw exposure each. We assume energy resolutions, energy ranges and effi- 
ciencies similar to the current versions of the experiments, and include backgrounds at target 
levels. Our analysis is based on full likelihood constructions for the experiments. We also 
take into account present uncertainties on hadronic matrix elements for neutralino-quark cou- 
plings, and on halo model parameters. We generate synthetic data based on four benchmark 
points and scan over the CMSSM parameter space using nested sampling. We construct 
both Bayesian posterior PDFs and frequentist profile likelihoods for the model parameters, 
as well as the mass and various cross-sections of the lightest neutralino. Future ton-scale ex- 
periments will help substantially in constraining supersymmetry, especially when results of 
experiments primarily targeting spin-dependent nuclear scattering are combined with those 
directed more toward spin-independent interactions. 

Keywords: dark matter theory, dark matter experiments, cosmology of theories beyond the 
SM, supersymmetry and cosmology 
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1 Introduction 

The Standard Model (SM) of particle physics (for an introduction, see ref. [1]) is superbly 
consistent with all existing laboratory experimental data. There are however various reasons 
to consider the SM as an effective theory valid only up to certain energies. These stem from 
conceptually irritating problems existent in the theoretical structure of the SM itself, and 
its inability to explain certain celestial observations. One of the most subtle issues in the 
first category is the so-called "gauge hierarchy problem". The problem is that the Higgs 
mass should be subject to very large (but experimentally excluded) quantum corrections, 
whose absence can only be explained in the framework of the SM by requiring extremely 
fine-tuned cancellations between physically-unrelated components of the theory. From the 
second category, the standard ACDM model of cosmology (for an introduction, see refs. [2, 
3]) is the strongest argument for new physics. The new physics called for ranges from a 
long-awaited reconciliation of Einstein's classical theory of gravity with the quantum-field- 
theoretical nature of the SM (required to describe very early moments of the Universe, and 
maybe also the properties of some extreme astrophysical objects) to the lack of any compelling 
explanation for dark matter (DM) or dark energy (see e.g. refs. [4-10]). If we add to this 
evidence the absence of any explanation for the origin of the gauge symmetries of the SM, 
the large differences in strengths of the gauge interactions, the number of matter generations, 
the origin of fermion masses, the strong CP problem, the specific scale and mechanism of 
electroweak symmetry breaking, and neutrino masses and mixings, going beyond the SM 
seems entirely unavoidable [1]. 
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Supersymmetry (SUSY), if realised at TeV energy scales, is arguably the most favoured 
theory beyond the SM. TeV-scale SUSY can solve many of the problems listed above, and 
pave the way for the resolution of many others in some broader theoretical framework (for 
an introduction, see refs. [11-13]). This weak-scale SUSY provides a natural solution to the 
hierarchy problem (by eliminating the quadratic divergences in the Higgs mass via the can- 
cellation of the contributions from SM particles and corresponding superpartners), suggests a 
natural mechanism for electroweak symmetry breaking (via renormalisation group evolution) 
and contains viable DM candidates (such as the lightest neutralino, gravitino, sneutrino and 
axino). Gauge-coupling unification can be achieved in models with SUSY [14] and many of 
the questions about the mathematical structure of the SM can be gracefully answered by as- 
suming some grand unified theory (GUT) above this unification scale, in which the SM gauge 
group is extended to some larger groups (e.g. SU(5) or 50(10)). Finally, SUSY, either as 
incumbent in the framework of supergravity (SUGRA) or as the essential ingredient of most 
versions of string theory, delivers extensive scope for unifying gravity and other fundamental 
interactions. 

Although it is not yet known how exactly SUSY is implemented in Nature (if at all), 
one can take a phenomenological approach and upgrade the SM to its SUSY extensions 
by simply adding new degrees of freedom and interactions to the SM Lagrangian. The 
most conventional framework for doing this is that of the so-called Minimal Supersymmetric 
Standard Model (MSSM), in which the particle content is increased by a minimum set of extra 
degrees of freedom necessary for supersymmetrising the model; i.e. every SM particle gets 
a "superpartner" . Were SUSY an exact symmetry of Nature, the structure of the MSSM 
would be highly restricted, including just one more parameter than the SM. This would 
however imply equal masses for SM particles and their superpartners, immediately excluding 
the theory by simple low-energy laboratory experiments. However, all the desirable features 
of SUSY remain entirely intact if it is broken around the TeV scales, in such a way that the 
gauge hierarchy remains stabilised, original gauge invariance of the model is preserved, and 
the renormalisability of the model is not violated. This is implemented in the MSSM by 
introducing the so-called "soft SUSY-breaking" terms to the (previously) SUSY-symmetric 
Lagrangian. 

The full MSSM Lagrangian contains more than a hundred free parameters to be deter- 
mined experimentally. The model is phenomenologically rich, but its huge parameter space 
makes any comparison of the model with real data highly challenging (see e.g. refs. [15-19] 
and references therein). There are several ways to approach this problem, so as to extract 
useful information about the structure of the model from experimental constraints. In gen- 
eral, one can pick out a particular SUSY-breaking mechanism (for some reviews, see e.g. 
refs. [20, 21]) that relates or even unifies many of the model parameters at certain ener- 
gies (usually much higher than the electroweak scale), and then study its phenomenological 
consequences at low energies by means of the renormalisation group equations (RGEs). Al- 
ternatively, one can impose experimentally-motivated relations directly on the parameters at 
low energy, so as to reduce their number. The latter proposal is justified by the fact that even 
small values of many of the soft parameters induce large flavour-changing neutral currents 
(FCNCs), or CP violation beyond that observed experimentally; these parameters can then 
be (more or less) reasonably set to zero. 

As far as SUSY-breaking mechanisms are concerned, several schemes have been proposed 
in the literature. These include gravity or Plank-scale mediation (PMSB) [22-27], gauge 
mediation (GMSB) [28-33] and extra-dimensional or anomaly mediation (AMSB) [34, 35]. 
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These schemes relate the model parameters in very different ways, giving rise to very differ- 
ent phenomenological predictions (for a comparison of some mediation models using existing 
data, see ref. [36]). A hybrid approach also exists, where one imposes phenomenological 
assumptions on model parameters at high energy scales, then studies low-energy predictions 
by solving the RGEs down to interesting scales. The Constrained MSSM (CMSSM) [37] 
is one of the most popular examples in this class of models, in which various boundary 
conditions are imposed at the GUT scale (^10 16 GeV). These conditions are motivated by 
minimal supergravity (mSUGRA) [22] as the simplest PMSB scenario. Although the number 
of free parameters is dramatically reduced by going from the MSSM to the CMSSM, many 
phenomenologically interesting features of SUSY remain intact, making the CMSSM partic- 
ularly interesting for phenomenologists trying to compare SUSY models with experimental 
data. 

SUSY signatures can be probed in a variety of complementary ways. The most obvious 
way is to explore physics at TeV scales directly at particle colliders; the Large Hadron Collider 
(LHC) will play a central role in this endeavour. If supersymmetric particles show up at such 
experiments, any positive measurement of their mass spectrum can put strong constraints 
on SUSY models and their parameters [38, 39]. Even if no positive signal is observed up 
to relatively high energies, one can still exclude extensive regions of many interesting SUSY 
scenarios using the resultant bounds on sparticle masses. 

There are however various other possibilities for testing SUSY. These are almost all 
related to the abilities of such models to provide DM candidates (for some general reviews, see 
e.g. refs. [40-43], and for an introduction to supersymmetric DM in particular, see ref. [44]). 
If DM consists of Weakly Interacting Massive Particles (WIMPs) thermally produced in 
the early Universe, these can be detected either directly, via interactions with SM particles 
in atomic nuclei of large-volume target materials on Earth, or indirectly by measuring the 
primary and secondary products of their self-annihilation in the sky. These products may take 
the form of photons, neutrinos, antimatter or other cosmic rays. To these detection channels 
we should also add the strong constraint on the relic abundance of DM, garnered from various 
ground-based telescopes and space missions. WMAP measurements of the cosmic microwave 
background anisotropics [4] have been perhaps the most influential of these. Direct and 
indirect searches are crucial for establishing that a particular candidate WIMP is in fact 
stable on cosmological timescales, and constitutes the majority of DM; whilst the LHC will 
be able to produce and discover heavy neutral particles, it will not be able to confirm their 
stability beyond the time they take to exit the detector [39]. 

Over the past few years, both direct and indirect detection experiments have seen huge 
boosts in both quality and quantity. In indirect detection (ID), gamma-ray telescopes and 
other particle detectors are currently observing potential products of DM annihilation from 
different astrophysical objects, through different annihilation channels. Some examples are: 
the Fermi gamma-ray space telescope [45], together with large air Cerenkov gamma-ray 
telescopes (ACTs) such as VERITAS [46], MAGIC [47] and H.E.S.S. [48] for detection with 
photons, the PAMELA satellite [49] together with Fermi, H.E.S.S. and some balloon missions 
such as ATIC [50] for detection with electrons, positrons and antiprotons, and IceCube [51] for 
detection with neutrinos. These experiments have already started touching upon interesting 
regions of SUSY parameter spaces (see e.g. refs. [18, 52]), but in some cases this ability 
should improve significantly in the next few years. Future experiments such as the AMS- 
02 space shuttle mission [53] and the GAPS long-duration balloon flight [54] provide the 
possibility of detecting antideuterons, another promising ID channel. The major problem 
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with ID searches is the large contamination of signals with backgrounds from astrophysical 
sources. This becomes especially problematic when one searches for a signal from the Galactic 
Centre. 

The situation is even more promising for direct detection (DD) experiments. Compared 
to the ID case, background contamination is much better understood. DD experiments are 
also relatively cheap and far easier to construct than ID ones. These factors have motivated 
several experimental groups to build such experiments. Currently, a copious number exist, 
sometimes with very similar properties, something that is always good for error reduction. 
These experiments have already reached sufficient sensitivity to probe small parts of SUSY 
parameter spaces, but this sensitivity is expected to significantly improve in the near future. 
It is therefore very interesting for both SUSY phenomenologists and DD experimentalists 
to get some feeling for how strong future DD limits and detections will be for SUSY DM 
models. It is our intention to provide that feeling in this paper. 

In this paper, we study the ability of future ton-scale DD experiments to constrain 
the parameter space of the CMSSM. For these purposes, we consider ton-scale versions of 
three current experiments (CDMS, XENON and COUPP), which are together expected to 
cover a large part of the remaining CMSSM parameter space. The analysis is performed 
using a state-of-the-art parameter estimation technique optimised for Bayesian statistical 
inference, and is based on full likelihood constructions for the experiments. We include typical 
energy resolutions, thresholds and efficiencies for such experiments, as well as backgrounds 
at target levels. We take into account different uncertainties in the values of Milky Way halo 
parameters, as well as the hadronic matrix elements for neutralino-nucleon couplings. 

The paper is organised as follows: in Sec. 2 we start with a brief review of the theory of 
direct detection of WIMPs, in particular the neutralino in the framework of SUSY DM. We 
describe our different assumptions regarding halo models and the calculation of neutralino- 
nucleon scattering cross-sections, as well as our treatment of different astrophysical and 
particle physics uncertainties. We then specialise the discussion to our chosen ton-scale 
experiments CDMS, XENON and COUPP, and present our likelihoods for the experiments. 
After discussing our benchmarks and generated data, we introduce our SUSY model and 
parameters, scanning strategy and statistical framework. Results are presented and discussed 
in Sec. 3, followed by concluding remarks and future prospects in Sec. 4. 



2 Framework, data and analysis 

2.1 Direct detection of supersymmetric WIMPs 

WIMP direct detection experiments seek to measure the energy deposited when a WIMP 
interacts with a nucleus in a detector [44, 55, 56]. If a WIMP of mass m scatters elastically 
from a nucleus of mass M, it will deposit a recoil energy E nr = {n 2 v 2 /M)(l — cos#), where 
ji = mMj (m + M) is the reduced mass of the WIMP-nucleus system, v is the speed of the 
WIMP relative to the nucleus, and 9 is the scattering angle in the centre of mass frame. The 
differential recoil rate per unit detector mass for a WIMP mass m, typically given in units 
of cpdkg -1 keV -1 (where cpd is counts per day), can be written as 

where q = \/2ME nr is the nucleus recoil momentum, a is the WIMP-nucleus cross-section 
in the limit of zero momentum transfer, F 2 (q) is a form factor due to the finite size of the 
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nucleus (normalised to F 2 (0) = 1), p is the local WIMP density, and information about the 
WIMP velocity distribution is encoded into the mean inverse speed r){E nr ,t), 



rj(E nr ,t) 



/(u,t) 



(2.2) 
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represents the minimum WIMP velocity that can result in a recoil energy E nr and /(u, t) is 
the (time-dependent) distribution of WIMP velocities u relative to the detector. 

Real experimental apparatuses cannot determine the event energies with perfect pre- 
cision. The recoiling nucleus (or recoiling electron, in the case of some backgrounds) will 
transfer its energy to either electrons, which may be observed as e.g. ionisation or scintilla- 
tion in the detector, or to other nuclei, producing phonons and heat. Measurements of these 
signals are then used to make an estimate E of the recoil energy for the event that produced 
those signals. However, due to intrinsic stochastic fluctuations in the chemical and physical 
processes that produce and propagate the ionisation, scintillation, and phonons, and due to 
technical limitations in measuring those signals, the measured signal (and thus estimated 
recoil energy E) exhibits some random variation for a given true recoil energy E nr , resulting 
in a finite (imperfect) energy resolution. In addition, limitations in detection efficiency and 
data cuts applied to reduce backgrounds mean that not all recoil events will be detected 
and/or accepted. Thus, in a real experiment, the expected observed spectrum is 



where (j){E, E nr ) AE is the probability of an event with recoil energy E nr being observed 
with a measured energy between E and E + AE (in the limit AE 0) after accounting for 
efficiencies (data cuts) and the finite energy resolution. There are two important distinctions 
between Eqn. (2.1) and Eqn. (2.4). (1) The latter is the observed rate, which is smaller than 
the former (nuclear recoil rate) as some events will fail to be observed and/or fall into the 
analysis window. (2) The latter equation is in terms of the measured energy E rather than 
the true recoil energy E nr and accounts for the random fluctuations in the measurements. In 
many cases, the fluctuations are small and taking E = E nr is a reasonable approximation. 
However, for many experiments, the relative size of the fluctuations grows as E nr becomes 
small, so the finite energy resolution should be taken into account whenever low energy recoils 
become important (typically the case for light WIMPs, as they can produce only low energy 
recoils). 

The expected number of observed (signal) recoil events between E m [ n and £ max is given 



where M^et is the detector mass, T is the exposure time, the sum is over the isotopes present 
in the detector material (possibly involving multiple elements), and fa is the mass fraction 
of each isotope. Here, we have neglected the time dependence in the source of which 
is described in Section 2.1.2. In principle, the time dependence can be properly accounted 
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for by replacing T in the equation above with an integral over the time that the experiment 
is running. However, the time dependence is small and not relevant for the experiments we 
consider in this paper, so we will use only the time averaged recoil spectrum of Eqn. (2.1). 

We use DarkSUSY [57] to calculate the differential recoil rate of Eqn. (2.1), albeit with 
some non-default parameters used in calculating the scattering cross-sections and describing 
the dark matter halo distribution. In the following sections, we examine these cross-sections 
and halo models in more detail. 

2.1.1 Cross-section 

For an MSSM neutralino, the only neutralino- quark couplings that do not vanish in the 
non-relativistic limit, and are thus not suppressed in the scattering of the low velocity relic 
neutralinos with baryons (the basis of direct detection), are the axial vector and scalar 
couplings [44, 58]: 

£> = «2gX7 /i 7 5 X^7/i7 5 ^ + asqXXQQ- (2-6) 

The axial vector and scalar couplings give rise to spin-dependent (SD) and spin-independent 
(SI) cross-sections for elastic scattering of a neutralino with a nucleus, respectively. Different 
form factors apply in each case, so that 

aF i {q) a SD Fi D (q) + a sl Fgi(q) (2.7) 

in Eqn. (2.1), where a SD and a SI are the SD and SI scattering cross-sections. The two 
cross-sections are described below. 

Spin-dependent cross-section (SD). The SD WIMP-nucleus cross-section can be writ- 
ten in terms of effective couplings to the proton (a p ) and neutron (a n ) as 

a s » = -G 2 F v 2 A 2 J(J+l), (2.8) 

7T 

where [i is the neutralino- nucleus reduced mass, J is the spin of the nucleus, 

A = -j (a p (S p ) + a n (S n )) , (2.9) 



and 1 



= v-^ A ( P ) an = V^A<"). (2.10) 



The factors Aq N ^ (N € {p, n}) parametrise the quark spin content of the nucleon and are 
only significant for the light (u,d,s) quarks. A combination of experimental and theoretical 
results tightly constrain the linear combinations [60] 



4 p) = A<p) - aJj P) = 1.2695 ± 0.0029 (2.11) 

and [61, 62] 



» = A (p) + a (p) _ 2A ( P ) = q 585 ± 0.025, (2.12) 



1 The effective spin-dependent couplings to the proton and neutron are sometimes instead given as and 
G™, differing from the above couplings only by the normalisation: Ga — 2\/2Gf cln [57, 59]. 
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while the strange quark contribution is given by the COMPASS result [63], 

A^ p) = -0.09 ± 0.01 (stat) ± 0.02 (syst) « -0.09 ± 0.03 , (2.13) 

where we have conservatively combined the statistical and systematic uncertainties. The 
proton and neutron spin contributions are related by an interchange of A u and (i.e. 
Al n) = A^, A^ n) = Al p) , and A^ n) = A^ p) ). The SD couplings a p and a n can be given in 
terms of the SD scattering cross-sections for a WIMP with the proton (o"p D ) and neutron 
(o"n D ), respectively, via Eqn. (2.8). 

The SD form factor F^ B (q) is dependent upon the spin structure within each nucleus and 
is thus unique to individual nuclei. We use the default form factors implemented in DarkSUSY, 
which are taken from Ref. [64] for iodine and xenon and from Ref. [65] for germanium (an 
interacting shell model form factor calculation is not implemented in DarkSUSY for fluorine, 
for which a simple exponential form factor is used). These and additional SD form factors 
may be found in the review of Ref. [66] . 

Spin-independent cross-section (SI). The SI WIMP-nucleus cross-section can likewise 
be written in terms of effective couplings to the proton (/ p ) and neutron (/ n ) as 

a SI = V[Z/ p + (A-Z)/ n ] 2 , (2.14) 

TV 

where fi is again the reduced mass and Z and A — Z are the number of protons and neutrons 
in the nucleus, respectively. For the scalar neutralino-quark couplings of Eqn. (2. 6), 2 

Jn_ = y 2 f (N) y as, ( } 

m N ^ JT « m a 27 JTG ^ m a ' y J 

g=u,d,s H g=c,b,t 

with 

mjv/g = (N\m q qq\N) (2.16) 

and 

g=u,d,s 

Here is the contribution of quark q to the composition of nucleus N. The contribution 

from the heavy quarks (c,b,t) arises through an anomaly involving gluons; the factor f^} is 
fixed by this anomaly [67, 68]. The matrix elements (p|^|p) for the light quarks (u,d,s) can 
be determined from estimates of the linear combinations [69, 70] 

X^v = i(m u + m d ) (p| (uu + dd) |p> = 64 ± 8 MeV (2.18) 



and [71-74] 
along with [70, 75] 



a = \(m u + m d ) (p| (uu + dd- 2ss) |p) = 36±7MeV, (2.19) 
(p| (dd - ss) |p) 



2 As with the spin-dependent case, the spin-independent couplings to the proton and neutron are sometimes 
given with an alternate normalisation: = 2/ p and = 2/ n [57, 59]. 
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Parameter Mean /s.d. 



Distribution Ref. 




-0.09 ±0.03 
36 ±7 
64±8 



normal 
log-normal 
log- normal 



[63] 
[71] 

[69, 70] 



po (GeV/cm 3 ) 



0.40 ±0.15 
235 ± 20 
235 ± 20 
550 ± 35 



log-normal 



[78-81] 



v rot (km/s) 
vo (km/s) 
v esc (km/s) 



normal 
normal 
normal 



[82-85] 



[86] 



Table 1. Likelihood distributions for parameters used to calculate the scattering cross-sections and to 
define the halo model. The parameters are defined in the text. The likelihoods for po and v ro t have been 
chosen based on a variety of estimates and do not represent any one particular estimate. The distribution in 
vo was chosen so that the allowed values of vo and v ro t cover a range of reasonable smooth dark matter halo 
models. Notably, we treat these two parameters as independent and do not make the common assumption 
that vo = v r ot. See the text for further discussion regarding the choice of distribution for each parameter. 

and the quark mass ratios ra u /m<i = 0.553 and m s /m<i = 18.9 [76]. The uncertainties in 
z and the quark mass ratios are negligible compared to the uncertainties in T, nN and ao- 
The pion-nucleon sigma term parametrises the nucleon mass shift from the chiral limit 
and is determined from pion-nucleon scattering, while the combinations of matrix elements 
corresponding to ao and z are determined from octet baryon mass differences. The neutron 
matrix elements are related to the proton matrix elements by an interchange of (7V|uu|7V) 
with (iV|dd|iV) (i.e. (n|uu|n) = (p|dd|p), (n|dd|n) = (p|uu|p), and (n|ss|n) = (p|ss|p)). 

As with the SD case, the SI couplings f p and f n can be given in terms of the SI 
scattering cross-sections for a WIMP with the proton (a^ 1 ) and neutron (o^ 1 ), respectively, 
via Eqn. (2.14). For an MSSM neutralino, Eqn. (2.15) typically yields f n ~ / p , so that 
a^ 1 ~ cTp 1 ; SI scattering is then typically characterised by the single parameter a^ 1 (aside from 
the neutralino mass). MSSM models often predict aff* » a^ 1 . However, the contributions 
to the total SI cross-section of a nucleus are a coherent sum over the individual proton and 
neutrons within, which yields the ~ A 2 scaling of the SI cross-section apparent in Eqn. (2.14). 
On the other hand, the contribution from the spins of individual nucleons tend to cancel in 
the SD case, so there is no such scaling with nucleus size for a SD . As a result of this 
scaling, heavier elements are often expected to have larger SI cross-sections than SD cross- 
sections, even if the reverse is true for a lone proton or neutron. For this reason, most direct 
detection experiments use heavy elements such as xenon, iodine, and germanium as their 
target materials and they mainly aim to observe or place constraints on SI interactions. 

DarkSUSY implements a variety of models for the SI form factor F^q). We use the 
default DarkSUSY form factors, which is usually the commonly used Helm form factor [56], 
but is a Fourier-Bessel or sum-of-Gaussians expansion for a limited number of isotopes in 
which these expansions are available. For the target elements used in our hypothetical ex- 
periments (see Section 2.2), most of the most significant carbon and germanium isotopes use 
Fourier-Bessel expansions (the exceptions being 13 C and 73 Ge), while the remainder use a 
Helm form factor. See Ref. [77] for a detailed discussion of these form factors. 

Parameters and likelihoods. We vary the quantities Ag P \ ao, and E^ in our scans with 
the likelihood distributions given in Table 1. For the positive parameter S^, we use a log- 
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normal distribution rather than a normal distribution as the latter would allow for unphysical 
negative values. In practice, though, there is little difference between the two distributions 
for the assumed mean and variance except in the tails of the distribution (beyond the 2-3a 
level) . We also avoid (likely disfavoured) small values of ao by using a log-normal distribution 
as these small values are suppressed in this distribution relative to a normal distribution. 3 As 
with S^jv, there would be little difference in practice had we used here a normal distribution 
for do as the two distributions are substantially similar in this case, except at their tails. 
We additionally require > ao, to avoid the strange scalar matrix element becoming 
(unphysically) negative. 

The parameters a^f\ and affect the determination of the SD scattering cross- 
section a SD for a given SUSY model. The uncertainties in these parameters as given in Table 1 
and in Section 2.1.1 result in variations in the calculated a SD by ~ xl.5, almost entirely due 
to the uncertainties in Ag P ^. As variations in and would not significantly affect our 
results, we fix them to their mean values and vary only Ag P ^. 

Uncertainties in both ao and yield significant variations in the predicted SI scat- 
tering cross-section a SI for a given SUSY model, by as much as an order of magnitude at 
the 2(j level. In many cases, the strange quark contribution in Eqn. (2.15) is the dominant 
contribution to Jn^ As (p|ss|p) oc T, nN — do, /n ~ (^ttn — ^o) 2 (approximately) in this 
case. The uncertainties in do and gain additional significance due to this dependence 
on the difference of the two terms, which has a larger total uncertainty than the individual 
parameters themselves. Given the very strong dependence of cr SI on the values of do and 

the uncertainties in these two parameters are the most important to take into account 
when fitting a SUSY model to a direct detection result. 

We neglect uncertainties in parameters other than these three (e.g. quark mass ratios) 
as the resulting uncertainties in the calculated SD and SI cross-sections are negligible in 
comparison. A more detailed discussion of how these and other parameters affect the elastic 
scattering cross-sections may be found in Ref . [87] . 

2.1.2 Halo model 

The dark matter halo in the local neighbourhood is likely to be composed mainly of a smooth, 
well mixed (virialised) component with an average density po ~ 0.4 GeV/cm 3 . The simplest 
model of this smooth component is the Standard Halo Model (SHM) [88]. The SHM is an 
isothermal sphere with an isotropic, Maxwellian velocity distribution characterised by an rms 
velocity dispersion a v . Sufficiently high velocity WIMPs would escape the Galaxy's potential 
well and will not be present in the halo in any significant quantities; the lack of these high 
velocity WIMPs can be accounted for by truncating the Maxwellian distribution at some 
escape velocity v esc , yielding 

fM = JN me v-», "~- v H for|v|< Uesc 

otherwise. 




To first order in SU(3) breaking, cfq ~ m r ^ yh (m« + ms — 2mjv), where rh = \ (m u + ma) [74, 75]; a 
negative value of ao would thus suggest the proton and neutron are not the lightest baryons. Even with higher 
order corrections, it would be difficult to rectify a negative or even small positive value of ao with observed 
baryon mass splittings. 

4 Due to the quark mass factor in Eqn. (2.16), fr s > fT u ,fr d - In the CMSSM and some other SUSY 
models, a^ q oc m q , so that the ^p- terms in Eqn. (2.15) are not suppressed for heavier quarks by the explicit 

quark mass appearing there. Combined with the ^ factor that suppresses the heavy quark contributions, the 
strange quark term will thus dominate unless a3 S is suppressed relative to other quark couplings. 
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Here 

7V esc = erf(» - 2zexp(-z 2 )/7r 1/2 , (2.22) 

with z = v esc /vo, is a normalisation factor that corrects for the truncation in the Maxwellian. 
Here we have also defined 

v = V2/3a v , (2.23) 

the most probable speed (in the non-truncated Maxwellian). 

The dark matter halo is essentially non-rotating, while the Sun moves with the disk, 
which rotates about the centre of the Galaxy at a speed v rot . The halo thus exhibits a bulk 
motion relative to Earth, so that 

f(u,t) = /(v obs (t) + u) (2.24) 

is the velocity distribution as seen from Earth's frame, where v b s (£) is the motion of the 
observer (Earth) relative to the rest frame of the dark matter halo described by Eqn. (2.21). 
The motion relative to this stationary halo is 

Vobs(^) = v L sr + v 05pec + V®(t) , (2.25) 

where vlsr = (0, ^ r ot 5 0) is the motion of the Local Standard of Rest in Galactic coordinates, 
v©, P ec = (10, 13,7) km/s is the Sun's peculiar velocity (see e.g. Refs. [89, 90] and references 
therein), and Ve(£) is the velocity of the Earth relative to the Sun. This last term varies 
throughout the year as the Earth orbits the Sun, leading to an annual modulation in the 
velocity distribution and, thus, the recoil rate of Eqn. (2.1) [88, 91]. However, v rot ~ 200- 
300 km/s, while the Earth's orbital speed is 30 km/s, so V© makes only a small contribution 
1° v b s . In addition, the Earth's orbital plane is inclined to the plane of the Galactic disk 
by 60° so that changes by ~ \V® (~15 km/s) due to the Earth's orbit. The result is 
that i; bs fluctuates by < 5% throughout the year. This results in a small time dependence 
of the velocity distribution that yields a small annual modulation in the recoil spectrum of 
Eqn. (2.1); observation of this modulation is the goal of experiments such as DAMA/Nal 
and DAMA/LIBRA [92, 93]. This small fluctuation is not important for the direct detection 
experiments and cases considered here; for simplicity, we drop the small V©(t) term in 
Eqn. (2.25) and use only the yearly average for v G b s . Further discussion of the annual 
modulation effect, along with analytical expressions for the integral in Eqn. (2.2) for this 
model, may be found in Ref. [94]. 

The isotropic, Maxwellian velocity distribution of Eqn. (2.21), intended to describe a 
class of smooth spherical halo models, should be considered only a first approximation of the 
local halo profile. Oblate/prolate or triaxial halos would be expected to have an anisotropic 
velocity distribution, though the assumption of isotropy is not necessarily a poor one if the 
level of oblateness/triaxiality is small, as is probably the case for the Milky Way. In addition, 
the formation of the Milky Way via merger events throughout its history leads to significant 
structure in both the spatial and velocity distribution of the dark matter halo. In the inner 
portions of the Galaxy, including the region containing the Sun, this structure is thought to 
have been substantially mixed to the point that a smooth distribution such as Eqn. (2.21) is a 
reasonable approximation. However, this does not preclude the possibility that a neutralino 
population from a late merger event (such as a tidal stream from a dwarf galaxy in the 
process of being absorbed by the Milky Way) might contribute an appreciable amount to the 
local density, with a velocity distribution significantly different from the well mixed portion. 
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However, the presence of such structure locally, while possible, is not probable and we ignore 
it here (see e.g. Ref. [95]). 

We note that the choice of excluding high velocity WIMPs by truncating the Maxwellian 
of Eqn. (2.21) at the local escape velocity is somewhat arbitrary. An alternative ad-hoc 
method of removing the high velocity WIMPs from the velocity distribution via subtraction 
instead of truncation can be seen in e.g. Refs. [96, 97]. Both of these forms are attempts to 
remove the most problematic part of the velocity profile "by hand." Chaudhury et al. [98] 
use King models to more self-consistently handle the finite size and mass of the Galaxy in 
determining the velocity distribution; in these models, the maximum WIMP velocity in the 
halo is, in fact, somewhat below the escape velocity. While a Maxwellian might reasonably 
approximate the bulk of the WIMP velocity distribution, the appropriate form of the tail 
is not so clear (but perhaps it will be clarified by numerical simulations [99]). Fortunately, 
the tail of the velocity distribution is unimportant for most direct detection experiments 
and WIMP candidates, the exception being very light WIMPs (^10 GeV or less), which 
we do not consider. 5 The imposition of a cut-off at the escape velocity v esc above serves to 
avoid unrealistic cases; however, any result that depends significantly on the value of this 
cut-off should be viewed with caution, as it likely also depends strongly on the (probably 
inaccurately) assumed form for the tail of the velocity distribution. 

Local dark matter density (po). The dark matter halo density has long been understood 
on the basis of Galactic modelling to have an average value of ^0.3 GeV/cm 3 in the solar 
neighbourhood, but with an uncertainty on the order of a factor of 2 [78]. This often-quoted 
value (though not the most recent or precise estimate) is typically taken as the fiducial case 
in determining direct detection constraints. As one can always write Eqn. (2.1) in terms of 
the quantity pq<j, direct detection experiments are unable to separately constrain the local 
density and scattering cross-section. Direct detection cross-section constraints are given 
either explicitly or implicitly in terms of the quantity ^a^ 1 (and similar for the SD scattering 
cross-sections), where £ = po/(0.3 GeV/cm 3 ). This allows for easy rescaling of constraints for 
other local density values, while the choice of a common reference value allows for comparison 
of constraints for different experiments. 

For scans over SUSY parameters, however, knowledge of the true local density is im- 
portant. If the assumed local density is too high or too low, the scans will find SUSY models 
yielding scattering cross-sections that are too low or too high, respectively. Thus, our limited 
knowledge of the true local density must be taken into account. 

Estimates of the local density are often made by fitting a variety of galactic observations 
to a particular model of the Galaxy that includes both dark matter and baryons. The 
reliability of these estimates are limited not only by statistical and systematic errors in the 
various observations, but in how well the assumed model matches the true Galaxy. Catena 
and Ullio [79] have suggested that for a spherical halo model, the local dark matter density 
can be very accurately determined as 0.385 ± 0.027 GeV/cm 3 for an Einasto density profile 
and 0.389 ± 0.025 GeV/cm 3 for a Navarro-Frank- White (NFW) profile. On the other hand, 
Weber and Boer, using a different method, find the local dark matter density is more weakly 
constrained, taking values of 0.2-0.4 GeV/cm 3 for spherical models and even higher for 
oblate models [80]. Salucci et al. have tried to estimate the local density using a more model 

5 In practice, there is typically very little difference in direct detection constraints when using the two 
different methods for modifying the Maxwellian distribution (truncation and subtraction). For CDMS or 
XENON exclusion curves, one finds the low mass fall-off of the constraint in the cross-section vs. mass plane 
is shifted by -1/2 GeV [100]. 
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independent method, finding p = 0.43 ± 0.11 ± 0.10 GeV/cm 3 [81]; see their paper for a 
description of the two uncertainties. Pato et al. have shown that, when relaxing complete 
spherical symmetry of the dark matter halo, the dark matter density near the stellar disk 
can be 1-41% larger than the spherically-shelled average at the same distance to the Galactic 
Centre [101]. This last issue should affect the estimates of Catena and Ullio, so the local halo 
density is probably not so precisely constrained as their stated values imply, though their 
novel method for determining the local density might also be applied to more complicated 
halo models. 

Velocity distribution. The velocity distribution characterised by Eqn. (2.21) is fixed 
by the three parameters v YOt , vq, and v esc . Based upon a variety of measurements, earlier 
estimates of the local disk circular velocity v rot placed it at ^220 km/s [82]. This 220 km/s has 
been the canonical value recommended by the International Astronomical Union (IAU) and 
is the most widely assumed value when determining direct detection constraints. 6 However, 
a recent estimate by Reid et al., based on Galactic masers, suggests a higher value of 254 ± 
16 km/s [83]. A reanalysis of the same data by McMillan and Binney find that v TO t ranges from 
200 ± 20 km/s to 279 ± 33 km/s, depending on how the Galactic rotation curve is modelled 
and what parameters are included in their fits, though the more extreme cases suggest a 
distance to the Galactic Centre i?o that is significantly lower or higher than typically found 
using other techniques [84]. The two most favoured models of McMillan and Binney yield 
235 ± 19 km/s and 258 ± 32 km/s, but a a major point of discussion is how estimates of the 
Sun's peculiar velocity v 05pec impact these fits. Bovy et al. have also reanalysed the same 
maser data to find v rot = 244 ± 13 km/s, or v rot = 236 ±11 km/s if other observations are 
included in the fit [85]. See Ref. [102] for an example of how different values of v rot impact 
direct detection constraints. 

The velocity dispersion cr Vl or most probable speed vq, of the dark matter halo is not 
directly observable, but it is related to the disk rotational velocity. For a spherical dark 
matter halo, vq ~ v YOt . For specific halo models, the relation is more concrete: for an 
isothermal spherical halo (the SHM), vo = v TOtl while for an NFW profile, vo would be ~10% 
smaller than v TO t [103] (the exact relation depends on the scales used in the NFW profile). 

WIMPs travelling faster than the local escape velocity for the Milky Way will escape 
the Galaxy and thus will not contribute to any direct detection signals. High velocity stars in 
the RAVE survey have been used to estimate the local escape velocity to be v esc = 544 km/s 
(498-608 km/s at the 90% confidence level), though this estimate is based upon a disk 
rotation speed of 220 km/s and could change for other values of v YOt [86]. 

Parameters and likelihoods. We assume a spherical halo model as described above 
and allow the nuisance parameters po 5 ^rot 5 and v esc to vary in our scans with likelihood 
distributions as given in Table 1. As noted previously, Eqn. (2.21) is only a first approxima- 
tion to the dark matter halo velocity profile; more general models of the halo would allow 
for triaxiality of the halo, anisotropy of the velocity profile, and the possible presence of 
structure remaining from recent merger events. We continue to use the spherical model for 
its simplicity, and note that our results will still be representative of how uncertainties in 
describing the smooth halo background affect the ability of scans to constrain the neutralino 
parameters. In any case, the direct detection experiments considered here are not likely to 

6 Often in the literature, the halo velocity distribution is characterised by a solar velocity of 232-233 km/s, 
which is obtained from this rotation speed of 220 km/s after including the Sun's peculiar velocity. See 
Eqn. (2.25) and the surrounding discussion. 
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be able to distinguish between more complicated models of the smooth halo and this first 
order approximation for SUSY models producing less than 0(100) events in the detectors. 7 

While many of these four parameters are highly correlated for a specific halo model 
(such as an isothermal sphere, where vq = ^ ro t) 5 we do not assume any particular halo model, 
so we take these quantities to be independent. This allows the scans to cover a wider range 
of realistic smooth halo backgrounds, allowing our lack of knowledge of the true halo profile 
to be taken into account in our reconstructions of CMSSM parameters. 

The distribution for the local density po (0.40 ± 0.15 GeV/cm 3 ) has been chosen to 
be representative of typical averages and uncertainties found in the literature, and does not 
correspond to any one estimate [79-81, 101]. A log-normal distribution is chosen to avoid 
inappropriately small local densities and to provide what is most likely to be an asymmetric 
distribution. The disk rotation velocity of v TO t = 235 ±20 km/s has been chosen to be similar 
to current estimates [83, 85], apart from the more extreme cases found in Ref. [84], and is 
still consistent with the older, IAU recommended estimate of 220 km/s [82]. For the velocity 
dispersion, we expect vq ~ ^ r ot 5 so the distribution in vo should have a relative variance 
at least as large as that of v rotl but the distribution should also allow for a ~10% relative 
variation between vq and v YOt that correspond to a variety of reasonable smooth halo models. 
We take vq = 235 ± 20 km/s, which satisfies both requirements. This choice for a central 
value induces a mild preference for an isothermal spherical halo (vq = ^ ro t) 5 but the variance 
still allows for velocity profiles consistent with reasonable alternative halo profiles (such as 
the NFW). Note that, even if the s.d. for v YOt had been much smaller, we would still wish to 
maintain the s.d. for vq of ^20 km/s in order to allow for a variety of halo models. For the 
escape velocity, we use v esc = 550 ± 35 km/s, comparable to the estimates made in Ref. [86]. 
We ignore the asymmetry in the distribution apparent in Figure 7 of Ref. [86] . We further 
note that this estimate was based on models of the Galaxy that assumed i; ro t = 220 km/s; we 
ignore the fact that this v^sc 

estimate should be adjusted for the different values of v TOt that 
we use. Again, we apply this distribution to indicate how uncertainties in halo parameters 
impact the scans and do not restrict ourselves to any particular halo model used to estimate 
the various halo parameters. In any case, the value of the escape velocity is not particularly 
important for the experiments we consider as signals in these experiments are dominated by 
contributions from neutralinos in the bulk of the velocity distribution rather than the tail 
(at least for the reasonably massive neutralinos we consider here). 

See Refs. [97, 103, 105] and references therein for further discussions regarding dark 
matter halo parameters and their impact on direct detection. See also Refs. [39, 106] for 
examinations of the impacts of these astrophysical uncertainties on reconstruction of WIMP 
parameters in direct detection results. 

2.2 Likelihoods for ton-scale CDMS, XENON and COUPP 

To examine the prospects for constraining SUSY models with future direct detection re- 
sults, we will use ton-scale extrapolations of the current experiments CDMS [107, 108], 
XENON10/100 [109-112], and COUPP [113, 114]; these (hypothetical) future detectors will 

7 The annual modulation signal and directional distribution of nuclear recoils searched for by experiments 
such as DAMA [92, 93] and DRIFT [104], respectively, are more sensitive to the halo velocity profile than 
experiments such as CDMS, XENON, and CoGeNT. These latter three experiments would be required to 
observe many events in order to distinguish the modest changes induced in Eqn. (2.2) by the use of a more 
complicated halo model. 
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Experiment 


Target 


Exposure 


Energy 


Event measurements: 






[kg-year] 


range [keV] 


Number? Energies? 


CDMS1T 


Ge 


1000 


10-100 


yes yes 


XEN0N1T 


Xe 


1000 


8-75 


yes yes 


COUPP1T 


CF 3 I 


1000 


>10 


yes no 



Table 2. Future ton-scale experiments considered in this work, based on extrapolations from current gener- 
ation experiments. The efficiencies and energy resolution of each experiment are discussed in the text. 

be referred to as CDMS1T, XENON1T, and COUPP1T, respectively. 8 In all three cases, 
we assume 1000 kg-years of raw exposure, with efficiencies, energy resolutions, and energy 
ranges similar to the present day versions of these experiments. The primary characteristics 
of the experiments are shown in Table 2. 

We define the likelihood of seeing N events with observed energies £i,..,jv m & direct 
detection experiment as 

C DD (N,E 1 _ N ) = P(N\vHot) J] f(Ei), (2.26) 

i=l. .N 

where P(N\fi tot ) is the probability of seeing N events for a Poisson distribution with average 
Utot and f(E) is the probability of an event having an observed energy E. The latter is 
normalised such that J^ max dE f(E) = 1, where ^min-^max is the energy range for the 
experiment. The quantity Utot and functional form of f(E) are dependent on the CMSSM 
parameters and cross-section nuisance parameters, whose dependences arise only through the 
neutralino mass and couplings to the nucleons (m, a p , a n , / p , / n ), as well as the halo nuisance 
parameters. For COUPP1T, which cannot determine event energies, only the number of 
events, we simply use the Poisson probability for the observed number of events 

C DD (N) = P(N\ntot) . (2.27) 

This special case will be discussed below. 

We briefly discuss here the use of Eqn. (2.26). This unbinned likelihood cannot be used 
to directly compare different experimental results (realisations) for a given set of CMSSM 
and nuisance parameters as the number of spectrum factors is iV, which is not constant. 
This means that this form has units of keV~ N and the likelihoods for two different experi- 
mental results (with different N) do not even have the same units. One way of comparing 
experimental realisations is to examine the likelihood ratio 

r= Cdd(N 1 E 1 _ n \0) 
C DD (N,E l9mm , N \0) ' 

where represents the SUSY and nuisance parameters and is the set of parameters that 
maximises the likelihood for that set of experimental results [115] (see Ref. [116] for an 
example of the usage of this ratio in analysing CDMS results). However, Eqn. (2.26) pro- 
vides all the dependence of the likelihood on the SUSY and nuisance parameters for a given 
experimental result and is sufficient for our purposes here. 

8 Ton-scale versions of these experiments are in various stages of planning, though the characteristics of 
those detectors may differ from what we use for our hypothetical detectors here. 



- 14 - 



(2.28) 



Cosmic rays, trace radioactive isotopes, and other sources of radiation induce recoils in 
a detector, providing a number of background events in addition to the WIMP-induced (sig- 
nal) recoil events. Many direct detection experiments aim to reduce the expected background 
to ^1 background event by applying a variety of techniques to both reduce the amount of 
radiation entering the detector, and to discriminate background from signal events. The for- 
mer typically involves reducing radioactive contaminants during detector fabrication, placing 
detectors deep underground to reduce the cosmic ray flux, placing shielding around detectors 
to reduce the influx of external radiation, etc. Discriminating background from signal is often 
done by observing two or more signatures from a recoil in the detector, such as ionisation and 
phonons (CDMS) or ionisation and scintillation (XENON). This allows very efficient removal 
of electron recoil events induced by /3 and 7 radiation, as recoiling electrons tend to produce 
the two signatures in a different ratio than recoiling nuclei do (where WIMP interactions in- 
duce only nuclear recoils, not electron recoils). Though this discrimination is efficient (often 
rejecting 99.9% or more of the 7-induced recoil events), the number of such background events 
can be in the 1000s, so some leakage into the signal region can still occur. Neutrons, on the 
other hand, induce nuclear recoils, which are indistinguishable from WIMP-induced nuclear 
recoils, and are harder to discriminate. However, the neutron radiation level is considerably 
smaller than 7-rays, and can be further reduced by removal of radioactive contaminants and 
using a muon veto to detect cosmic ray showers (which can produce neutrons). In addition, 
neutrons can scatter multiple times in a detector; observation of a multiple scatter is a clear 
signature of a neutron, because the probability of a WIMP scattering more than once is 
negligible. This allows all multiple-scattering events to be excluded from further analysis. 

Though the goal of reaching a background level of ^1 event might seem overly ambitious, 
it is not an unreasonable one. The current generation of detectors have been able to reach this 
target, such as CDMSII [108] (2 background/signal events observed), XENON100 [112] (0 
events), and COUPP [114] (3 events). Though extending this background level to much larger 
detectors is no easy task, we note that the discrimination efficiency tends to increase with the 
size of the detector. The larger detector volume increases the likelihood of seeing neutrons 
multiple-scattering within the detector and thus allows for better rejection of neutrons. A 
larger volume also allows for self-shielding: external radiation is unlikely to penetrate far 
into the active volume without interacting; ignoring recoil events in the outer portion of the 
detector volume (referred to as a fiducial volume cut) can significantly reduce the background 
level with only a modest reduction in the expected WIMP signal, as WIMPs easily penetrate 
the detector volume. 

In addition to the WIMP-induced recoil spectrum, we assume a background spectrum 
consisting of a component that is flat in energy and a component that decreases exponentially 
with increasing energy 

^ = C f + C e e- E '*> , (2.29) 

where Cf and C e are normalisation constants chosen such that the expected number of back- 
ground events from the flat and exponential components are hbj and /iB je , respectively. The 
above form was chosen to reflect the two common energy dependences of typical backgrounds: 
some are fairly independent of energy (flat spectrum), while some backgrounds tend to in- 
crease rapidly at low (and decreasing) energy. The particular form of the background for any 
given experiment depends on the detection technique and analysis cuts; some experiments 
may expect to have backgrounds more consistent with only one of the components, while 
some may have a background spectrum not adequately approximated by Eqn. (2.29) (e.g. 
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there may be a Gaussian peak due to a particular radioactive decay). We take the above 
form for all the experiments we consider, as a representative example of what a background 
spectrum might be, not as an indication of what the actual spectrum will be if/when such 
experiments are built. Indeed, while the expected number of background events hb can often 
be found for direct detection experiments, the spectra of these backgrounds can rarely be 
found in the literature even for existing experiments. 9 We presume the future experiments 
reach their low background level goal and take hbj — I^B,e — 1-0 in all three cases, with Eq 
= 10 keV. The total expected number of (signal+background) events is then 

l^tot = us + , (2.30) 

where [is is the expected WIMP signal of Eqn. (2.5) and \ib = Hbj + VB,e is the total 
expected background. The expected (signal+background) energy spectrum is 

dRtot _ dRs_ dR B , v 

dE ~ dE + dE ' [Z ' 6i) 

where is the WIMP spectrum of Eqn. (2.4). With this spectrum, the spectral factor in 
Eqn. (2.26) is 

1 dR tot [ E ™*„ dR tot , OQO s 

f( E ) = Jf,dE-' Wlth Nf = J E . dE ^- (2 " 32) 

We now describe the CDMS1T, XENON1T, and COUPP1T experiments we use in our 
study. We emphasise that, while their characteristics are loosely based on their existing 
present-day versions, we make various assumptions and approximations that do not neces- 
sarily reflect the true behaviour of the current experiments or their future incarnations. Our 
intent is to define experiments that provide a variety of targets and realistically account for 
such issues as finite energy resolution. 

2.2.1 CDMS 1-ton 

The Cryogenic Dark Matter Search (CDMS) experiments use germanium detectors to observe 
both phonons and ionisation produced in recoils [107, 108]. For our hypothetical future ton- 
scale CDMS-like experiment (CDMS1T), we take a 1000 kg- year germanium exposure over 
energies of 10-100 keV. 

The 10 keV threshold is typical of the CDMS low background analyses. Though CDMS 
is capable of examining events at energies much lower than 10 keV, background discrimination 
rapidly becomes worse at lower energies and the spectrum at these energies is significantly 
contaminated with backgrounds [117, 118]. The 10 keV threshold already allows for signifi- 
cant sensitivity to most WIMPs; inclusion of lower energy events is generally only important 
when searching for very light WIMPs (less than ^10 GeV) incapable of producing any events 
above threshold. Due to the increased background contamination, addition of the low-energy 
spectrum to the analysis would give only a modest or even negligible improvement to the 
sensitivity to WIMPs heavier than ~ 20 GeV. As we do not focus on low mass WIMPs, 
the low-energy portion of the spectrum is not important for our analysis and this 10 keV 
threshold is sufficient for our purposes. Due to the form factor and mean inverse velocity 



9 The lack of background spectra may be the result of large systematic uncertainties that significantly affect 
its determination. Nevertheless, if /when a potential WIMP signal is detected in an experiment, characterising 
the background contributions will be of utmost importance in order to use the signal to constrain SUSY 
parameter spaces. 
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factor 77 in Eqn. (2.1), the recoil rate at energies of ~ 100 keV and above is far lower than 
near threshold, so even though CDMS is capable of detecting such recoils, these high energy 
events contribute very little to the overall rate (at most 5%). An upper limit on energies 
is imposed rather than including this small high-energy signal, as backgrounds can easily 
produce a similar or greater number of events at these energies. 

We now define the factor in Eqn. (2.4) for CDMS1T. The CDMS detectors allow for 
a very good energy resolution; we use [119] 

a 2 (E nr ) = (0.293) 2 + (0.056) 2 (£ nr /keV) [keV 2 ] (2.33) 

to describe the normally distributed fluctuations of the observed energy E about the true 
recoil energy E nr . This energy resolution technically applies to low energy electron recoil 
events, though we assume here that nuclear recoil events have a similar resolution. Though 
the efficiency e of the various cuts is in principal energy dependent, we assume a constant 
30% efficiency, or 10 

£(£ nr ) = 0.3, (2.34) 

which is comparable to the 25-32% efficiency over 10-100 keV in the latest CDMS analysis 
[108]. With this energy resolution and efficiency, 

<KE, E nr ) = e{E nr ) 1 = e -(E-E nr) y2aHE nr ) (2>35) 
^2i\a l [E nr ) 

2.2.2 XENON 1-ton 

The XENON10 [109-111] and XENON100 [112] experiments use a dual phase liquid/gas 
xenon detector and observe both scintillation and ionisation from recoil events. For our 
hypothetical future ton-scale XENON-like experiment (XENON1T), we take a 1000 kg- year 
xenon exposure over energies of 8-75 keV. The XENON experiments have a poorer energy 
resolution than CDMS, but the heavier target (A « 131 for Xe versus A « 73 for Ge) 
gives a somewhat better sensitivity to SI scattering due to the A 2 scaling. The same issues 
as with CDMS, regarding the energy range, also apply here. Lower energy events can be 
examined, but to reach much lower energies requires a different type of analysis, leading to 
greater background contamination [120]. As with CDMS, ignoring lower energy events only 
significantly impacts the XENON sensitivity to light WIMPs and is not important for our 
purposes. 

As with CDMS, we define a normally distributed fluctuation of E about E nr , with [121] 



a(E nr ) = (0.579 keV) y/E nr /te*V + 0.021£ nr , (2.36) 
and a constant efficiency 

e(E nr ) = 0.3, (2.37) 

which is the approximate efficiency for XENON100 when including the fiducial volume cut 
[112]. With this energy resolution and efficiency, the factor <j) for XENON1T is given by 
Eqn. (2.35). This energy resolution is a fairly simplified form for the XENON10 experiment 
and the resolution at low recoil energies is actually much more complicated for XENON- 
like experiments [122]. However, the above form still describes an experiment with poorer 



10 Depending on how the cuts are defined, the efficiency may alternatively be defined as a function of the 
observed energy E, i.e. s(E), rather than as a function of the true nuclear recoil energy E nr . 
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energy resolution than CDMS, and should introduce the same type of relative limitations 
in constraining parameters as would be expected from a XENON-like experiment. We note 
that there are currently significant issues in determining the energy calibration of XENON- 
like experiments (i.e. estimating the energy E from a given scintillation signal), which affect 
both the threshold and energy resolution determinations; see Refs. [123, 124] and references 
therein for further discussion. We assume that these energy calibration issues will be resolved 
by the time a ton-scale xenon experiment is built, and neglect any systematic errors that 
would otherwise be introduced in the analysis of our XENON1T experiment. 

2.2.3 COUPP 1-ton 

The Chicagoland Observatory for Underground Particle Physics (COUPP) is a bubble cham- 
ber experiment, different from CDMS and XENON in that it does not measure event energies, 
only the number of events above a threshold nuclear recoil energy [113, 114]. The target ma- 
terial, CF3I, provides some sensitivity to SI scattering, primarily through the heavy iodine 
nuclei, but the high number of nuclei with non-zero spins (all the fluorine and iodine nu- 
clei) allows for a significant sensitivity to SD scattering. The sensitivity to SD scattering is 
likely to be much higher than in CDMS or XENON, allowing COUPP to potentially con- 
strain WIMP parameters in a manner different than (and complementary to) the other two 
experiments. 

The principle behind the operation of the COUPP bubble chamber begins with its use 
of superheated liquid CF3I within the chamber. When sufficient energy is deposited in a 
small region, a bubble forms in this liquid, undergoing runaway growth and causing the 
liquid to boil. However, if insufficient energy is deposited or if the energy is spread out over 
a larger region, it is thermodynamically favourable for a bubble to recollapse rather than 
expand, and the liquid remains in its superheated state (such bubbles are too small to be 
seen). Electron recoils, induced by (3 or 7 radiation, deposit energy over too large a region 
to induce runaway bubble growth, and are thus not observed. Nuclear recoils, on the other 
hand, deposit their energy over a small region and, if above a certain minimum recoil energy, 
will generate an observable bubble in the detector. This minimum energy (threshold) is 
adjustable by changing the operating temperature and pressure of the liquid. Essentially 
all nuclear recoil events above threshold produce bubbles, which are indistinguishable and 
give no indication of their recoil energies (other than that they exceeded the threshold). 
Thus, COUPP provides no direct measurement of the recoil energy spectrum. The number 
of observed bubbles is the only observation that can be used to constrain WIMP properties, 
leading to the use of Eqn. (2.27) instead of Eqn. (2.26) for the likelihood. Alpha decays 
unfortunately also lead to bubble formation, and constitute a source of background events. 
These alpha decays produce a different acoustic signature in the detector than nuclear recoils, 
allowing for some background discrimination; see Ref. [114] for details. 

For our hypothetical ton-scale COUPP-like experiment (COUPP1T), we take a 1000 kg- 
year CF3I exposure with a threshold of 10 keV. We take an efficiency of 



comparable to the efficiency of the most recent COUPP run after fiducial volume and acoustic 
signature cuts [114]. As energies are not measured, there is no finite energy resolution to 
take into account, and we take E — E nr by setting 



e(E nr ) = 0.7 



(2.38) 



0(E,E nr ) =e(E nr )S(E-E nr ). 



(2.39) 
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Carrying out the integration in Eqn. (2.4) yields a spectrum ^ proportional to the true 
nuclear recoil spectrum of Eqn. (2.1), though ^ will only be used in Eqn. (2.5) and does 
not represent the "observed" spectrum here, as no energy spectrum is actually observed. 

In principle, some sensitivity to the energy spectrum can be gained by dividing the 
COUPP exposure into multiple periods with different thresholds. The relative rates observed 
at different thresholds is dependent upon the recoil spectrum; however, sensitivity to the 
spectrum is far weaker via this technique than with experiments like CDMS and XENON. 
In practice, each period at a different threshold would be treated as a separate "experiment" 
with a likelihood defined by Eqn. (2.27). We do not consider this case here. 

2.3 Benchmarks and generated data 

In order to assess the ability of our hypothetical ton-scale experiments CDMS1T, XENON1T 
and COUPP1T to constrain the CMSSM, we first look at Fig. 1, where the strongest limits 
on WIMP properties from current direct detection (DD) experiments are given. The black 
and red curves, shown in the m-o-a^ 7 plane, are exclusion limits at the 90% confidence level 
from the CDMS and XENON experiments under the assumption of a standard local halo 
configuration. We see that almost all WIMPs with a^ 1 > 10 -7 pb are excluded. 

So, what will our hypothetical ton-scale DD experiments add to our knowledge of 
WIMPs in the CMSSM? One way to examine this is to assume that such experiments do 
not see any signal events (i.e. no WIMPs exist with cross-sections within the reach of these 
experiments), and then put stronger exclusion limits on WIMP properties. For these types of 
analyses one can use the likelihoods we constructed for different experiments, and scan over 
the model parameter space by setting the observed number of signal events to zero. This will 
exclude regions of parameter space that predict some signal events. What we are interested 
in however is how well the experiments can zero-in on parts of the parameter space if they 
see an actual signal. For this purpose, we need to choose some hypothetical 'benchmark' 
models, and test how well we can recover them in our hypothetical experiments. 

We perform this investigation by selecting some benchmark points in the m-o-a^ 7 plane 
of Fig. 1 in such a way that they represent WIMPs with interesting properties that are not 
yet ruled out. For each point we randomly generate some hypothetical DD data based on the 
likelihoods of Sec. 2.2, and use that data in CMSSM scans to constrain the model parameter 
space. Our selected benchmark points are shown in Fig. 1 with red circles. Values of the 
neutralino mass, SI and SD scattering cross-sections for our chosen points are explicitly listed 
in Table 3. 

We have chosen four benchmarks (shown as BM1, BM2, BM3 and BM4 in Fig. 1). 
With benchmarks 1, 3 and 4 we cover the most interesting regions in the m-o-cr^ 7 plane with 
relatively high cross-sections. The points correspond to WIMPs with low, intermediate and 
high masses (~ 100 GeV, - 500 GeV and - 1100 GeV, respectively). These are WIMPs 
for which we expect a significant number of signal events in ton-scale experiments. We also 
select a point at low mass (~ 100 GeV) and relatively small cross-section (~4x 10 -11 pb), 
for which we expect only a couple of events above background. The cross-sections for this 
benchmark point approximately correspond to the sensitivity limits of our detectors; below 
this cross-section we do not expect them to be able to distinguish WIMPs from backgrounds. 
We do not look at higher masses with low cross-sections, because in those cases the expected 
number of events is so low that the detectors would see only backgrounds. For each of the 
points we have selected in Figure 1, we have chosen for the corresponding benchmark an 
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Figure 1. Strongest experimental limits on the spin-independent scattering cross-section of the neutralino 
and a proton cr^ 1 , as a function of the neutralino mass, from CDMS and XENON experiments. All WIMPs 
above the curves are excluded at the 90% confidence level if a standard local halo configuration is assumed. 
Red circles show the benchmark points chosen in our analysis. For corresponding benchmark values of spin- 
dependent cross-sections, see Table 3. 



arbitrary set of CMSSM parameters that yields the desired neutralino mass and SI cross- 
section, as shown in Table 3, in order to generate reasonable SD cross-sections. 

The neutralino mass and cross-sections for each of our benchmarks were determined 
from a particular set of CMSSM parameters. However, the same mass and cross-sections in 
each case can be approximately reproduced with many different CMSSM models, in some 
cases with very different parameter values, due to various degeneracies between and/or weak 
dependencies on some of the CMSSM parameters. Thus, direct detection experiments will 
only be able to directly constrain the neutralino mass and cross-sections and will be unable to 
significantly constrain all of the CMSSM parameters (m^ is the exception as it is strongly 
related to the neutralino mass). We therefore focus our studies on how direct detection exper- 
iments constrain the neutralino mass and cross-sections rather than on individual CMSSM 
parameters, though we do include the constraints for some of the latter. 

Table 4 shows the total and fractional numbers of expected events for the benchmark 
points, as well as the actual numbers we found in our randomly-generated data, based on 
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Quantity 


BM1 


BM2 


BM3 


BM4 


af (pb) 
(Pb) 
(Pb) 

m^o (GeV) 


3.9xl(T 9 
2.8xl(T 5 
2.1xl0- 5 
108.32 


3.9xl0- n 
1.5xl(T 8 
2.8xl(T 8 
95.51 


1.7xl(T 9 
4.8xl(T 7 
3.9xl0- 7 
513.45 


4.1xl0- 9 
4.5xl(T 7 
3.4xl(T 7 
1104.05 


m (GeV) 
m 1/2 (GeV) 
A (GeV) 
tan /3 


3354.7 
271 
3515.3 
40.2 


1773.8 
239.7 
4340.2 
22.3 


2008.8 
1194.7 
2840.6 
39.5 


3905.5 
2490 
3609.5 
22 



Table 3. (Top) Values of spin-independent (SI) and spin-dependent (SD) cross-sections of the neutralino 
as well as the neutralino mass at our benchmark points BM1, BM2, BM3 and BM4. (Bottom) CMSSM 
parameters for the benchmark points. These parameters were arbitrarily chosen to produce the desired 
neutralino mass and SI cross-section for each benchmark. We note that there are other choices for CMSSM 
parameters, sometimes with very different values, that can also produce the same neutralino mass and cross- 
sections in each case. 



the likelihoods. Our synthetic data set also includes randomly-generated recoil energies for 
the events, according to the expected energy spectrum. 

Given some scattering cross-sections and a neutralino mass, the expected and generated 
number of events, and their corresponding recoil energies, also depend on the assumed values 
for the halo parameters po, ^rot, ^0 and v esc . In generating our synthetic data, we set these 
values to the means of the probability distribution functions given in Sec. 2.1.2; values are 
given in Table 5. 

2.4 CMSSM scans 

2.4.1 Model, parameters and ranges 

The model of choice in our analysis is the CMSSM [37]. 11 Here one unifies all scalar, gaugino 
and trilinear coupling parameters of the MSSM to mo, m^ 2 and Ao, respectively. This 
unification happens at the gauge coupling unification (or GUT) scale (^10 16 GeV). Assuming 
such boundary conditions, one obtains low-energy quantities such as the masses of the SM 
particles and their superpartners by solving the RGEs from the GUT scale down to scales 
where these quantities are defined. One key assumption in the CMSSM is that electroweak 
symmetry is broken at TeV scales; this requirement fixes the magnitude of the Higgs/higgsino 
mass parameter \i in the MSSM superpotential. The sign of \i remains as a free parameter 
to be determined by experiment. In addition to mo, m^ 2 , A§ and sgn/x, the CMSSM 
possesses one more free parameter: tan/3, the ratio of the up-type and down-type Higgs 
vacuum expectation values. Compared to the SM, the CMSSM has therefore five new free 
parameters in total, four continuous and one discrete. Throughout this paper we only work 
with positive fi. We use SOFTSUSY [126] (available from ref. [127]) to solve the RGEs for 
every set of CMSSM parameters. 

There are also some SM quantities whose values have not been fully determined by 
experiment. The ones that impact the CMSSM predictions most are m^, the top quark mass, 
777,5, the bottom quark mass, a e m 5 the electromagnetic coupling constant and a 5 , the strong 

11 An alternative analysis using a purely phenomenological setup has been performed in ref. [125]. Results 
of that analysis are in generally good agreement with those we present in this paper. 
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BM1 


BM2 


BM3 


BM4 


CDMS1T 


l^tot 


135.2 


3.4 


18.7 


21 


Ps 


133.2 


1.4 


16.7 


19 


..SI 

Ps 


131.8 


1.4 


16.7 


19 


SD 
Ps 


1.39 


2.2xl0~ 3 


6.9xl0~ 3 


2.9xl0~ 3 


PB 


2 


2 


2 


2 


Nobs 


139 


4 


18 


23 


XENON1T 


Vtot 


204.1 


4.1 


24.6 


27.3 


Ps 


202.1 


2.1 


22.6 


25.3 




195.3 


2.1 


22.6 


25.3 




6.7 


l.lxlO- 2 


3.2xl0~ 2 


1.3xl0- 2 




2 


2 


2 


2 


N obs 


213 


4 


24 


28 


COUPP1T 


l^tot 


389.9 


5 


34.4 


38 


Us 


387.9 


3 


32.4 


36 




268.6 


2.9 


31.9 


35.7 




119.3 


6.6xl0~ 2 


0.56 


0.25 


Hb 


2 


2 


2 


2 


Nobs 


392 


5 


35 


38 



Table 4. Expected and generated numbers of events for our hypothetical ton-scale direct detection ex- 
periments CDMS1T, XENON1T and COUPP1T at the four benchmark points BM1, BM2, BM3 and BM4 
employed in our scans of the CMSSM parameter space, ^tot, Ms, /xf 7 , /i| D , \jlb and N Q bs denote the to- 
tal expected numbers of (signal+background) events, expected numbers of signal events, expected numbers 
of spin-independent signal events, expected numbers of spin-dependent signal events, expected numbers of 
background events and the total generated numbers of events in our particular realisations, respectively. 



Po 


^rot 


^0 ^esc 


(GeV/cm 3 ) 


(km/s) 


(km/s) (km/s) 


0.4 


235 


235 550 



Table 5. Values of halo parameters used for generating direct detection data at our benchmarks; these are 
the same as the mean values for distributions given in Sec. 2.1.2. 



coupling constant. We therefore add them to the set of free parameters in our scans. Here 
we evaluate the bottom quark mass m^ at m^ and denote it as mb(rrib) MS while for the 
top quark mass we use its pole mass and show it simply as m t . The electromagnetic and 
strong coupling constant s ar e evaluated at the Z-boson pole mass mz and are denoted as 
a ein (mz) MS and a s (mz) MS , respectively. MS means that these quantities are computed in 
the modifie d m inimal su btra ction renor mali sation scheme. Values of the four SM quantities 
mt, mb(mb) MS Oi ein (mz) MS and a s (mz) MS are strongly constrained by electroweak precision 
measurements and we therefore include them in our scans as nuisance parameters. We 
model the likelihoods for these parameters by normal distributions with means and standard 
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Observable 


Mean value 


Standard deviation 


Ref. 


m t (GeV) 


172.6 


1.4 


[128] 


m b (m b ) MS (GeV) 


4.20 


0.07 


[60] 


a s (m z ) MS 


0.1176 


0.002 


[60] 


l/a em (mz) MS 


127.955 


0.03 


[129] 



Table 6. Experimental constraints on SM nuisance parameters used in our analysis. 



Parameter 


Lower limit 


Upper limit 


m (GeV) 


60 


4000 


m 1/2 (GeV) 


60 


4000 


A (GeV) 


-7000 


7000 


tan f3 


2 


65 


m t (GeV) 


163.7 


178.1 


m h {m h f^ (GeV) 


3.92 


4.48 


a s (m z ) MS 


0.1096 


0.1256 


l/a em (mz) MS 


127.846 


127.990 




-0.18 





a (MeV) 


20 


65 


S^(MeV) 


40 


100 


p (GeV/cm 3 ) 


0.1 


1.2 


v rot (km/s) 


175 


295 


v (km/s) 


175 


295 


v esc (km/s) 


445 


655 



Table 7. Ranges of model and nuisance parameters employed in our scans. 



deviations given in Table 6. 

Given a point in the eight-dimensional parameter space of the CMSSM+SM nuisance 
parameters, one can theoretically predict properties of the neutralino, such as its mass and 
scattering cross-sections. However, as pointed out in Sec. 2.1.1, the cross-sections depend 
also on the values we assume for the neutralino-quark couplings through the hadronic matrix 
elements. As we discussed before, the uncertainties on three of these quantities, Ag P \ ao and 
Xt^/v, can significantly impact predictions for the cross-sections. We therefore vary them in 
our scans as cross-section nuisance parameters with likelihoods given in Sec. 2.1.1. 

Even though the neutralino mass and cross-sections can be calculated just from CMSSM 
parameters, SM nuisances and cross-section nuisances, the DD data observed by a detector 
depend also on the halo model and parameters (Eqns. 2.1 and 2.2). As discussed in Sec. 2.1.2, 
we do not fix the values of our four halo parameters po 5 ^rot 5 and v esc in our scans; they 
form yet another set of nuisance parameters with likelihoods given in Sec. 2.1.2. 

With four CMSSM parameters, four SM nuisances, three cross-section nuisances and 
four halo nuisances, the parameter space that we scan over becomes fifteen-dimensional. 
Ranges for all parameters used in our scans are listed in Table 7. 

The full likelihood for our model is a composite likelihood comprised of the likelihoods 
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for our ton-scale DD experiments and the likelihoods for the SM, cross-section and halo 
nuisance parameters. In addition, we require some physicality conditions at each sample point 
in the parameter space: (i) self-consistent solutions to the RGEs must exist, (ii) electroweak 
symmetry breaking (EWSB) conditions must be fulfilled, (hi) all particle masses must be 
non-tachyonic, and (iv) the neutralino must be the lightest supersymmetric particle (LSP). 
To ensure that all these conditions are satisfied, we disregard any points that do not fulfil 
these criteria simply by assigning them extremely small likelihoods. 

2.4.2 Statistical frameworks and scanning technique 

In order to make any inference about the parameters of a theoretical model by comparison 
with experimental data, one has to first decide which statistical framework to use, Bayesian or 
frequentist statistics. SUSY parameter estimation is no exception; for a general introduction, 
see e.g. ref. [130], and for some applications in SUSY phenomenology, see refs. [16, 17, 19] 
and references therein. 

In Bayesian statistics (for its applications in physics, see e.g. ref. [131], and for re- 
views of its applications in cosmology, see e.g. refs. [132-135]), probabilities are defined in 
a 'subjective' way and consequently, one can assign probabilities to model parameters. The 
cornerstone of Bayesian inference is Bayes' Theorem through which our 'prior' degree of be- 
lief about different values of model parameters is updated when new experimental data are 
taken into account. Our updated knowledge, given in terms of a multi-dimensional probability 
density function (PDF) and called the 'posterior' PDF, is then used to make statistical state- 
ments about the model and its parameters, including the construction of intervals/regions 
in the parameter space with certain levels of confidence ('credible' intervals/regions). This 
construction can be easily performed by integrating (or 'marginalising') the full posterior 
PDF over all unwanted parameters resulting in joint 'marginal' PDFs for the parameters of 
interest. 

In frequentist statistics, on the other hand, probability has an 'objective' definition, so 
probabilities cannot be assigned to model parameters. Instead of posterior PDFs, frequen- 
tists work only with the likelihood and the idea of repeated experiments. One widely-used 
method for constructing 'confidence' intervals/regions in the parameter space in the frequen- 
tist framework is the 'profile likelihood' construction [136, and references therein]. In this 
method, in order to make inference about some parameters, we maximise (or 'profile') the 
full multi-dimensional likelihood function along the unwanted parameters. For complex and 
poorly-understood models like the CMSSM, with non-Gaussian likelihoods, the profile like- 
lihood method is only an approximation to the exact frequentist construction of confidence 
regions proposed by Neyman [137]. One method that provides exact confidence regions is the 
so-called 'confidence belt' construction [115]; this method is however harder to implement 
numerically. 

Results of Bayesian and frequentist statistics, in particular the credible or confidence 
intervals/regions obtained, do not agree in general. One central ingredient of Bayesian in- 
ference is the set of prior assumptions one employs in the inference. This can significantly 
affect the final results, in particular when data are not constraining enough and/or the model 
at hand is too complex (for a detailed discussion of the impacts of priors for SUSY param- 
eter estimation, see ref. [16]). Bayesian statistics therefore provides a good measure of the 
robustness of a fit; if the fit depends strongly on priors, it is not to be considered definitive. 
On the other hand, results of frequentist inference including the profile likelihood framework 
are in principle independent of priors. This is however the case only if we have an accurate 
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construction of the likelihood function. How true this is depends on how powerful the em- 
ployed sampling algorithm is. Most existing scanning algorithms are optimised for Bayesian 
inference (for an example of frequentist scanning techniques, see e.g. ref. [17]), so any map- 
ping of the likelihood function based on these scans can be strongly impacted by prior effects 
(for detailed discussions, see refs. [17, 19]). 

The other interesting distinction between marginal posteriors and profile likelihoods is 
that the marginalisation procedure makes the posterior sensitive only to the total posterior 
probability mass in the subspace that has been marginalised over. Consequently, fine-tuned 
regions with spike-like likelihood surfaces cannot be probed. Profile likelihoods, however, are 
strongly sensitive to fine-tuned regions by construction, as confidence regions are constructed 
as iso-likelihood contours around the best-fit point. Both marginal posteriors and profile 
likelihoods should be considered in any scan, so as to acquire a complete picture of the 
favoured parameter regions; this is our strategy in this paper. 

In order to map the posterior PDFs and profile likelihoods for our SUSY model, we 
employ nested sampling [138, 139] to explore the parameter space. We use the MultiNest 
algorithm [140, 141], which is widely used in SUSY parameter estimation. MultiNest is 
one scanning option in SuperBayeS [16, 142, 143] (available from ref. [144]), a package for 
constraining CMSSM parameters using different types of experimental data. Nested sampling 
is designed and optimised for computing the Bayesian evidence (and hence, posterior PDFs) 
in a relatively fast and efficient way. One other popular class of scanning algorithms are 
Markov-Chain Monte Carlos (MCMCs), which are also optimised for Bayesian inference. 
Results of nested sampling and MCMC scans are very similar in the CMSSM [16], although 
the former requires remarkably less computational power. 

3 Results and discussion 

3.1 Bayesian credible and frequentist confidence regions 

Before we present the results of our scans when DD data are employed, we give in Fig. 2 
the two-dimensional marginal posteriors and profile likelihoods when only the priors and 
physicality constraints are taken into account. Plots are in the SI and SD scattering cross- 
sections of the neutralino a^ 7 , a^ D and a^ D versus the neutralino mass m^o, as well as the 
CMSSM parameter planes mo-m^ 2 and Ao-tan/3. We also show the same plots for the halo 

and hadronic mass matrix nuisance parameters po, v TO t, vo, ^esc, Ag P \ ao and T, nN . These 
show how our scanning algorithm samples the parameter space when no experimental data 
are provided, so can be thought of as the "effective priors" we place on the parameters by 
our choice of physicality constraints and scanning algorithm. The inner and outer contours 
in each panel (in dark and light blue for marginal posteriors and in yellow and red for 
profile likelihoods) represent 68.3% (la) and 95.4% (2a) confidence levels, respectively. Black 
dots show the posterior means while black crosses show best-fit (highest likelihood) points, 
respectively. 

These types of plots are always necessary when one wants to examine how constraining 
certain sets of data (from DD experiments in our case) are. A combination of both posterior 
PDFs and profile likelihoods show that in our scans some regions of the parameter space 
are sampled better than others. This becomes clear when we examine the plots for the 
DD quantities cr^ 7 , a^ D ', a^ D and m^o in particular. Obviously, the algorithm tends to 
sample low cross-sections and low neutralino masses at a higher resolution. For a detailed 
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Figure 2. Two-dimensional marginalised posterior PDFs and profile likelihoods for the mass and nuclear 
scattering cross-sections of the lightest neutralino, as well as the halo and cross-section nuisance parameters, 
when no direct detection data are included in the scans. The same plots are also shown for CMSSM parameters 
(last row). These plots essentially show the effects of explicit priors imposed on the parameters, as well as 
the physicality constraints. The inner and outer contours in each panel represent 68.3% (la) and 95.4% 
(2a) confidence levels, respectively. Black dots and crosses show the posterior means and best-fit points, 
respectively. Benchmarks are marked with green stars. 
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Figure 3. Two-dimensional marginalised posterior PDFs for the mass and nuclear scattering cross-sections 
of the lightest neutralino for benchmark 1 when different combinations of direct detection likelihoods are used 
in our scans. The inner and outer contours in each panel represent 68.3% (la) and 95.4% (2a) confidence 
levels, respectively. Black dots and crosses show the posterior means and best-fit points, respectively, and 
benchmark values are marked with green stars. 



discussion of these sampling effects, their relationship to the choice of priors and implications 
for statistical inference, see the companion paper [19]. For the halo and cross-section nuisance 
parameters, the plots are basically reflections of what we assumed for their likelihoods, i.e. 
normal and log-normal distributions. 

We then present our main results when synthetic data from the DD experiments CDMS IT, 
XENON1T and COUPP1T are included. They are shown in Figs 3-10 for different bench- 
marks. For each benchmark, we first present two-dimensional marginal posteriors for the 
three scattering cross-sections a^ 7 , a^ D and a^ D versus the neutralino mass m^o. Again, as 
in Fig. 2, la and 2a contours are shown in dark and light blue, respectively. Corresponding 
two-dimensional profile likelihoods are given in separate figures for the same planes of ob- 
servables. Here la and 2a contours are depicted in yellow and red, and posterior means and 
best-fit points are denoted by black dots and crosses, respectively. In all plots, we also mark 
the locations of the benchmark values by green stars. 

In each figure, we show the credible/confidence regions for different combinations of 
likelihoods for DD experiments: the first column in the left corresponds to the case where only 
CDMS1T is used, the second is for XENON1T only, the third corresponds to a combination of 



-27- 





Profile likelihood 




Flat priors 




CMSSM, |i>0 




CDMS only 



Profile likelihood 
Flat priors 
CMSSM, ji>0 
XENON only 



Profile likelihood 
Flat priors 
CMSSM, |i>0 
CDMS+XENON 



Profile likelihood 
Flat priors 
CMSSM, |i>0 
CDMS+XENON+COUPP 



mo (TeV) 



m^o (TeV) 



m x o (TeV) 



mo (TeV) 







Profile likelihood 






Flat priors 






CMSSM, \i>0 




w 


CDMS only 


( ' 

\1 





0.5 1 1.5 

m o (TeV) 



Profile likelihood 
Flat priors 
CMSSM, n>o' 
XENON only 



0.5 1 1.5 

m o (TeV) 



Profile likelihood 
Flat priors 
CMSSM, n>0 
CDMS+XENON 



0.5 1 1.5 

m o (TeV) 



Profile likelihood 
Flat priors 
CMSSM, n>0 
CDMS+XENON+COUPP 



0.5 1 1.5 

m o (TeV) 







Profile likelihood 






Flat priors 






CMSSM, |i>0 






CDMS only 







0.5 1 1.5 

mo (TeV) 



Profile likelihood 
Flat priors 
CMSSM, ji>0 
XENON only 



0.5 1 1.5 

mo (TeV) 



Profile likelihood 
Flat priors 
CMSSM, |i>0 
CDMS+XENON 



0.5 1 1.5 

mo (TeV) 



Profile likelihood 
Flat priors 
CMSSM, |i>0 
CDMS+XENON+COUPP 



0.5 1 1.5 

mo (TeV) 



Figure 4. As in Fig. 3, but for two-dimensional profile likelihoods. 



CDMS1T and XENON1T, and the last column shows the regions when all three experiments, 
CDMS1T, XENON1T and COUPP1T, are taken into account. Regions for COUPP1T alone 
were not determined as, without a measurement of event energies (the recoil spectrum), 
COUPP1T has almost no ability to constrain the neutralino mass and provides much weaker 
overall constraints than CDMS IT or XENON1T alone. However, as will be discussed below, 
COUPP1T can prove to be very useful in constraining various quantities when combined 
with one or both of the other experiments. 

In Fig. 11, we also present two-dimensional marginal posteriors and profile likelihoods 
for the CMSSM parameters mo and Aq versus m^ 2 and tan/3, respectively, when different 
benchmarks are used in our reconstruction procedure. Since the focus of this paper is the 
dark matter mass and cross-sections, we show only the cases where combination of all three 
DD likelihoods is used in the scans so as to avoid the large number of figures in the paper (for 
plots with other combinations of DD likelihoods, see the first preprint of the paper available 
on the arXiv). 

In Figs. 3 and 4 for benchmark 1, with relatively high cross-section and low neutralino 
mass, all the experiments constrain the combination of a^ 1 and m^o quite successfully, even 
when uncertainties on the halo and cross-section parameters are taken into account. Both 
the best-fit point and the posterior mean are located very close to the benchmark point in all 
cases. This is what we would expect from all three experiments when there are many signal 
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Figure 5. As in Fig. 3, but for benchmark 2. 



events (as in this case), most of which are SI events. 

By looking at the a p D -m^o and a^ D -m^o planes, we see that the capability of CDMS1T 
and XENON1T to constrain SD cross-sections is much less than for the SI case. This is to 
be expected, as these experiments do not have optimal target nuclei for SD detection. The 
situation is different for COUPP1T (rightmost column of Figs. 3 and 4). SD events for 
COUPP1T constitute a large fraction of the total events, which strongly breaks the degen- 
eracy in the a^-m-o and a^-m-o planes. This can be seen in the substantial narrowing 
of contours when data from COUPP1T is included in the fit. To clarify this picture even 
further, we show in Figs. 12 and 13 two-dimensional posterior PDFs and profile likelihoods 
in terms of the two effective couplings of the neutralino to the proton and neutron, a p and 
a n (Eqn. 2.10). These two quantities are strongly correlated in the CMSSM. By looking 
at the panels in the first row we see that for the cases where results from COUPP1T are 
not included, virtually no constraint is placed upon a p and a n , and \a p \ and \a n \ can take 
values as small as zero. With COUPP1T, however, relatively large values of \a p \ and \a n \ 
are favoured, meaning that in this case, SD interactions are important and COUPP1T has 
been able to detect them. Here, CDMS1T and XENON1T are able to constrain both the 
neutralino mass and the SI cross-section, but are unable to place any significant limits on the 
SD cross-sections due to the statistically insignificant number of SD events expected (see Ta- 
ble 4). Once the mass and SI cross-section are approximately determined, which COUPP1T 
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Figure 6. As in Fig. 4, but for benchmark 2. 



cannot do alone, a clear excess of COUPP1T events, relative to the expected number of SI 
events, becomes apparent: 392 events are observed in this pseudo-experiment, but only ^270 
SI events are expected for the neutralino mass and SI cross-section determined by CDMS IT 
and/or XENON1T. This ^120 event excess can be attributed to SD scattering, allowing the 
SD cross-sections to be constrained. In this case, COUPP1T proves to be a quite useful ad- 
dition to the other experiments and this situation illustrates how multiple DD experiments 
can provide complementary (rather than redundant) results. 

The prior effects seen in Fig. 2, and discussed previously, are strongly compensated for 
in this case and the credible/confidence regions are dominated by data rather than priors 
(either real or effective). The mo-m^ 2 plots (first column in the left of Fig. 11) show that 
DD data for benchmark 1 with low neutralino mass can strongly constrain m^; this is not 
surprising, given the strong correlation between the neutralino mass and m^2 in the CMSSM. 
For the Ao-tan/3 planes, DD data exclude regions at very high or very low Aq. Such data do 
not put strong constraints upon tan/3, except for the case where COUPP1T is added to the 
combination. Adding COUPP1T also excludes regions with very high or very low tan/3 (as 
shown in Fig. 11). 

For benchmark 2, Figs. 5 and 6 and Fig. 11 (second column in the left) show that the 
credible/confidence regions are strongly dominated by prior effects. For example, comparing 
panels in the first row of Fig. 5 for a^ 1 versus m~o with the corresponding plot in Fig. 2 
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Figure 7. As in Figs. 3 and 5, but for benchmark 3. 



shows that including data has only excluded high cross-sections at low neutralino masses. 
The situation improves slightly when the COUPP1T likelihood is added to the analysis; in 
this case, parts of the region at low cross-sections and high masses are also excluded. The 
reason why our DD experiments cannot constrain the model well in this case is the small 
number of signal events compared to background (Table 4), giving rather low statistics. 
The same conclusions can be drawn from the mo-m^ 2 and Ao-tan/3 plots for benchmark 
2 in Fig. 11. They are again dominated by priors except for very low values of Aq in the 
Ao-tan/3 plane, which seem to be excluded. Regions with low values of both mo and m^ 2 , 
as well as high-m^ regions are also somewhat disfavoured. The reason for the unusual 
a n -a p regions for the CDMSlT-only case with benchmark 2, shown in Figs. 12 and 13, is not 
completely clear, but we suspect it arises due to this particular pseudo-experimental result, 
which had four events at energies of 13.3, 79.7, 90.5, and 95.3 keV. The three higher energy 
events are likely due to a ~ 2a random upward fluctuation in the number of flat background 
events (an average of one flat background event is expected). These few high-energy events, 
with a lack of low-energy events, are not particularly consistent with either a neutralino 
spectrum or the exponentially-falling background. The combination of the few events and 
unusual energy distribution may have driven the scans to prefer atypical sets of parameters. 
We note that, though these regions are unusual, the true parameters still lie within the la 
contours. 
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Figure 8. As in Figs. 4 and 6, but for benchmark 3. 



Figs. 7 and 8 and Fig. 11 (third column in the left) present our results for benchmark 3. 
These plots show that compared to benchmark 1, which has a much lower neutralino mass 
(Figs. 3 and 4 and second column of Fig. 11), the experiments are less constraining in this 
case. Combining different experiments does not improve the results significantly. This is 
because the number of total events is largely reduced and most of the events come from SI 
interactions (Table 4). The latter point explains why the degeneracy in the two SD couplings 
is not broken in this case when the COUPP1T likelihood is taken into account. 

Finally, similar plots are given in Figs. 9 and 10 and Fig. 11 (last column) for benchmark 
4. The neutralino mass for this point is relatively large (~ 1 TeV). Our credible/confidence 
regions in this case are very similar to the previous case for benchmark 3. 

One interesting feature of Figs. 9 and 10 is that contrary to the other benchmarks, the 
true values for the mass and cross-sections do not fall within the la and 2a regions. From 
a Bayesian point of view, such a situation can naturally happen in the presence of strong 
priors that dominate the statistical inference, in the absence of high statistics provided by 
data. In a frequentist framework, however, where the inference should be independent of 
priors, any confidence regions and intervals provided by a successful parameter estimation 
procedure, 'must' show "proper coverage". The latter simply means that in an ideal case, 
regions with 95.4% (2a) confidence level, for instance, contain the "true" parameters in 95.4% 
of the times that a large number of scans are performed using synthetic data generated from 
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Figure 9. As in Figs. 3, 5 and 7, but for benchmark 4. 



a psudo-experiment. This indicates that the plots in the first row of Fig. 10, where is 
plotted against m~o for the two-dimensional profile likelihoods, and the benchmark is outside 
the 2a regions, can be perfectly acceptable as long as such plots do not appear in more than 
4.6% of the times. This is not the case however, as is discussed in detail in the companion 
paper [19], where statistical coverage is studied in SUSY parameter estimation when data 
from DD experiments are used (also see ref. [145] for a similar coverage study of the CMSSM, 
but in a different context). By looking at the first rows of Fig. 10, we observe that contours 
taper off at high masses, and only a few points exist in those regions even though very good 
fits continue to be found as neutralino masses increase. For neutralinos much heavier than 
the nuclear targets, as is the case here, DD experiments are unable to provide an upper 
limit on the neutralino mass and the best-fit contours in the cr^-m-o plane should extend 
indefinitely along a line of constant a^/m^o at high masses. This then suggests that our 
CMSSM scans, due to some "sampling effects", map out only the lower mass end of these 
expected contours. As discussed in ref. [19], we believe that this under-coverage seen in the 
profile likelihood results, is due to sampling effects caused by strong prior effects that reduce 
the ability of our scanning algorithm (i.e. MultiNest) to explore regions at high masses and 
cross-sections. MultiNest is a technique optimised for Bayesian statistics and, consequently, 
any frequentist inference it provides can be indirectly (and in some cases strongly) influenced 
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Figure 10. As in Figs. 4, 6 and 8, but for benchmark 4. 



by the prior effects. The configuration we use for our MultiNest scans (i.e. the particular 
values of the tolerance parameter and the number of live points) is the same as what is used 
in, for example, [16]. A few scans with a smaller value of the tolerance parameter and a 
higher number of live points (as suggested in ref. [146]) show that one might obtain better 
coverage in some cases (for more details, see ref. [19]). 

In order to see what the predictions of our scans are for other types of experiments, 
in Figs. 14 and 15 we give two-dimensional marginal posteriors and profile likelihoods for 
the velocity-averaged neutralino self-annihilation cross-section (cry), as a function of the 
neutralino mass m^o. These plots are commonly used when indirect detection constraints 
are put upon WIMPs like the lightest neutralino. Every row in the figure corresponds to one 
benchmark point when different combinations of the DD experiments CDMS1T, XENON1T 
and COUPP1T are employed. For panels in the first row of each figure, where benchmark 
1 (with relatively large scattering cross-sections and low neutralino mass) is studied, we see 
that our scans predict rather tight contours, especially when COUPP1T is considered. The 
neutralino mass is highly constrained, as we have seen in previous figures (Figs. 3 and 4). In 
addition, the neutralino annihilation and scattering cross-sections are strongly correlated, so 
because the latter are highly constrained for benchmark 1, the correlation results in relatively 
tight constraints on the annihilation cross-section. We suspect that a large fraction of the 
favoured region in this case is the 'funnel' region, where one of the neutral Higgs particles in 



-34- 





Figure 11. Two-dimensional marginalised posterior PDFs and profile likelihoods for CMSSM parameters 
and for all four benchmarks, when the combination of all three direct detection likelihoods is used in the 
scans. The inner and outer contours in each panel represent 68.3% (la) and 95.4% (2a) confidence levels, 
respectively. Black dots and crosses show the posterior means and best-fit points, respectively, and benchmark 
values are marked with green stars. 



the MSSM has roughly twice the mass of the lightest neutralino, and neutralino annihilation 
is increased by a mass resonance with that Higgs particle (for a discussion of distinct regions 
in the CMSSM parameter space and their relevance for indirect detection analysis, see e.g. 
ref. [18]). The funnel region seems to be found also in our scans for benchmark 2 (regions at 
low neutralino masses in second rows of Figs. 14 and 15). This funnel region is often difficult 
to find when employing the same priors on sparticle masses as we have here (flat priors), so 
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Figure 12. Two-dimensional marginalised posterior PDFs for the effective spin-dependent neutralino- 
neutron (a n ) and neutralino-proton (a p ) couplings, when different combinations of direct detection likelihoods 
are used in our scans. The inner and outer contours in each panel represent 68.3% (la) and 95.4% (2a) con- 
fidence levels, respectively. Black dots and crosses show the posterior means and best-fit points, respectively. 
Benchmarks are marked with green stars. 



it is interesting that they show up in our scans. 
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Figure 13. As in Fig. 12, but for two-dimensional profile likelihoods. 



We finally show in Table 8 experimental sensitivities (detection prospects) for the differ- 
ent benchmarks. The table gives the fraction of the time that our ton-scale DD experiments 
are expected to detect a dark matter signal at the 3 or 5a level based only on the number 
of events observed. We consider WIMPs to be detected if the observed number of events 
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Figure 14. Two-dimensional marginalised posterior PDFs for the velocity- averaged neutralino self- 
annihilation cross-section (av) as a function of the neutralino mass rn^o, when different combinations of 
direct detection likelihoods are used in our scans. The inner and outer contours in each panel represent 68.3% 
(la) and 95.4% (2a) confidence levels, respectively. Black dots and crosses show the posterior means and 
best-fit points, respectively. 



greatly exceeds that expected due to backgrounds alone. This happens when the observed 
number has a p- value smaller than 2.7 x 10 -3 and 5.7 x 10 -7 (corresponding to the 3a and 5a 
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Figure 15. As in Fig. 14, but for two-dimensional profile likelihoods. 



detection, respectively) for a Poisson distribution based on an average number of background 
events /ig. The average background for the individual experiments is /ig = 2.0 and this cor- 
responds to observing 8 (13) or more events at the 3a (5a) level detection. We also give in 
the last row the combined results of all three experiments. This provides greater statistical 
power than each experiment alone. For an average background \ib — 6.0 in this case, a 3a 
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Experiment 


BM1 


BM2 


BM3 


BM4 


CDMSIT 

XENONIT 

COUPPIT 

(Combined) 


1.000 (1.000) 
1.000 (1.000) 
1.000 (1.000) 
1.000 (1.000) 


0.023 (0.00006) 
0.057 (0.00034) 
0.133 (0.00202) 
0.275 (0.00940) 


0.998 (0.931) 
1.000 (0.996) 
1.000 (1.000) 
1.000 (1.000) 


1.000 (0.975) 
1.000 (0.999) 
1.000 (1.000) 
1.000 (1.000) 



Table 8. Fraction of the time that the direct detection experiments are expected to detect a dark matter 
signal at the 3a (5a) level for each benchmark model based only on the number of observed events. Dark 
matter is considered to be detected at the 3a (5a) level if the observed number of events greatly exceeds 
that expected due to backgrounds alone, specifically if the observed number of events has a p- value smaller 
than 2.7 x 10 -3 (5.7 x 10 -7 ) for a Poisson distribution based on an average number of background events hb- 
For the average background \ib — 2.0 assumed for each of the individual experiments, this corresponds to 
observing 8 (13) or more events. The last row represents the combined results of all three experiments, which 
provides greater statistical power than each experiment alone. A 3a (5a) detection in this case corresponds 
to observing 15 (22) or more total events for an average background \ib — 6.0. The use of event energies in 
addition to the number of events could give a moderate improvement in the detection prospects. 

(5(J ) detection corresponds to observing 15 (22) or more total events. Although we do not in- 
clude event energies in our appraisal of the detection capability of the experiments, this could 
however moderately improve the detection prospects. Our results show that benchmarks 1, 
3 and 4 will be detected in almost all cases while for benchmark 2 the best case would be 
when all three experiments are combined, with 27.5% and 0.94% detection probabilities at 
the 3 or 5a level, respectively. 

3.2 Effects of halo and cross-section nuisances 

In our analysis, we included uncertainties on the halo parameters po 5 ^rot 5 and v esc , as 
well as on the neutralino-nucleon coupling parameters Ag P \ ctq and E^, by treating them 
as nuisance parameters. It is interesting to know how large these nuisance effects are, and 
how strongly they impact our scanning results. In Figs. 16 and 17 we show the results of 
some scans for benchmark 3 with data from all three experiments, CDMSIT, XENONIT 
and COUPPIT, and where none or only some of these nuisances are included (for brevity 
purposes, we again only show the results in terms of the mass and scattering cross-sections of 
the neutralino, and refer the reader to the first preprint of the paper on the arXiv for similar 
plots in the CMSSM planes). 

By comparing panels in the first rows (the a^-m-o planes) we see that when the 
nuisance parameters are allowed to vary, the favoured region becomes wider in cross-section. 
This effect is mostly caused by the halo nuisances (compare the third and fourth panels in 
the first rows), and fixing those whilst only varying the cross-section parameters (second 
panels) does not change the size and shape of the confidence/credible region significantly, 
compared to the case where all nuisances are fixed (fourth panels). The latter behaviour is 
expected: changing the values of ao and E^ do not affect the value of a^ 1 that best fits the 
experimental results, but it does affect which CMSSM models yield that SI cross-section. An 
interesting feature in the case where all nuisance parameters are fixed is the sharp cut-off 
in the SD cross-sections a^ D and a^ D at high values (fourth panels in the second and third 
rows of Figs. 16 and 17); this can be understood from the correlation between SD and SI 
cross-sections. 
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(a) CS/H Nuis: 



(b) CS Nuis Only: 



(c) H Nuis Only: 



(d) No CS/H Nuis: 
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Figure 16. Two-dimensional marginalised posterior PDFs for the mass and nuclear scattering cross-sections 
of the lightest neutralino, for benchmark 3 when all, some or none of the halo and cross-section nuisance 
parameters vary along with the CMSSM parameters in our scans. CS and H in the subcaptions stand for 
cross-section and halo nuisances, respectively. Again, the inner and outer contours in each panel represent 
68.3% (la) and 95.4% (2cr) confidence levels, respectively. Black dots and crosses show the posterior means 
and best-fit points, respectively, and benchmark values are marked with green stars. 



An interesting case to look at is when our DD experiments highly constrain the CMSSM 
parameter space (i.e. benchmark 1), and when the likelihoods from all three experiments are 
used. Fig. 18 compares the results for two cases: (i) when all the halo and cross-section 
nuisances vary and (ii) when all the halo and cross-section nuisances are fixed to their mean 
values. This figure clearly shows that including COUPP1T in the latter case allows us 
to reconstruct the true values of the scattering cross-sections and neutralino mass almost 
perfectly, with uncertainties coming essentially only from our ignorance about the nuisance 
parameters. Fig. 18 also shows how well one can constrain the CMSSM parameters in two 
cases (i) and (ii). Although these measurements impose strong constraints on the values 
of m^ 2 , still wide ranges of values for mo, Aq and tan/3 are equally consistent with the 
data. This therefore indicates that other types of experiments, such as indirect detection or 
accelerator searches are required to break these degeneracies and pin down the actual values 
of the model parameters. 
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Figure 17. As in Fig. 16, but for two-dimensional profile likelihoods. 



Throughout the paper, in generating synthetic DD data at each benchmark point, we 
fixed the halo and cross-section nuisance parameters at their mean values given in Table 1. 
We believe that our results would not change significantly if we chose other benchmark values 
for these parameters within e.g. their 1 or 2a ranges and generated the simulated data using 
those new benchmark values of nuisances; this is essentially the whole idea behind marginal- 
isation. To demonstrate this, we scanned over the CMSSM parameter space assuming a 
set of randomly generated values for the nuisances: A^ p) = -0.124, a = 33.0 (MeV), 

= 57.6 (MeV), p = 0.308 (GeV/cm 3 ), v TOt = 203.6 (km/s), v = 238.7 (km/s), 
^esc — 538.2 (km/s). In other words, for each benchmark point (with CMSSM parame- 
ters given in Table 3), instead of the mean values of Table 1, we used the above values for 
nuisance parameters, and then generated the synthetic data to be used in the scans. Likeli- 
hood distributions used in these scans had the same forms (with the same mean values) as 
before (i.e. introduced in Table 1). We show the results in Fig. 19, where the two-dimensional 
posterior PDFs and profile likelihoods are compared to the case where the nuisances are fixed 
to their mean values. As expected, credible/confidence regions are not affected strongly; they 
are only slightly larger in the case where the true values of the nuisances differ from the means 
of their distributions. 
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(a) CS/H Nuis: 
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Figure 18. Two-dimensional marginalised posterior PDFs and profile likelihoods for the mass and nuclear 
scattering cross-sections of the lightest neutralino, as well as the CMSSM parameters, for benchmark 1 when 
all or none of the halo and cross-section nuisance parameters vary along with the CMSSM parameters in our 
scans. CS and H in the subcaptions stand for cross-section and halo nuisances, respectively. The difference 
between this figure and Figs. 16 and 17 is the choice of benchmark point; here, we consider a case where direct 
detection data strongly constrain the parameters (benchmark 1), whereas in Figs. 16 and 17 we considered a 
case where constraints were weaker (benchmark 3). Again, the inner and outer contours in each panel represent 
68.3% (la) and 95.4% (2a) confidence levels, respectively. Black dots and crosses show the posterior means 
and best-fit points, respectively, and benchmark values are marked with green stars. 
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(a) Mean values: 



(b) Random values: 



(c) Mean values: 



(d) Random values: 
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Figure 19. Comparison of two-dimensional marginalised posterior PDFs and profile likelihoods for two 
different choices of the 'true' halo and cross-section nuisance parameters. All plots are for benchmark 3. First 
and third columns show the results when the nuisances are set to their mean values when generating synthetic 
data. Second and fourth columns show the same, but for randomly chosen values of the nuisances. 



4 Summary and conclusions 

We have studied the ability of future ton-scale dark matter direct detection experiments to 
put constraints upon the parameter space of the CMSSM. The experiments we considered 
were ton-scale extrapolations of CDMS, XENON10/100 and COUPP. We assumed 1000 kg- 
years of raw exposure for all three experiments, with efficiencies, energy resolutions and 
energy thresholds similar to their current versions. In our construction of the likelihoods, 
we included the number of events seen by the experiments, as well as the observed event 
energies for CDMS 1-ton and XENON 1-ton. We modelled different background events with 
a background spectrum containing a flat and an exponentially-falling component. 

We generated synthetic data for some benchmark points with particular nuclear-scattering 
cross-sections and neutralino masses. We scanned over the CMSSM parameter space and re- 
constructed the marginalised posterior PDFs and profile likelihoods. These gave us Bayesian 
credible and frequentist confidence regions for different combinations of parameters and ob- 
servables. We have also taken into account uncertainties in different halo model parameters 
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and neutralino-quark couplings, by considering them as nuisance parameters and marginal- 
ising over them. 

Our results show that by combining all three experiments, one can successfully deter- 
mine the WIMP mass and cross-sections for significantly low-mass and high-cross-section 
neutralinos. This ability lessens substantially when WIMPs have much higher masses and 
lower cross-sections. We have shown that COUPP can be very helpful for breaking degenera- 
cies between different parameters when a large fraction of its signal event rate comes from SD 
interactions. We also observed that uncertainties in the halo model have a reasonably large 
effect on the constructed credible/confidence regions, whereas cross-section uncertainties do 
not change the results significantly. 

Finally, our results imply that even in the best-case scenario where the neutralino has 
relatively low mass and high cross-sections, direct detection experiments would not constrain 
all CMSSM parameters strongly; the tightest constraint would be on the gaugino mass pa- 
rameter, which is strongly correlated with the neutralino mass. This indicates that if indeed 
SUSY is realised in Nature, other types of experiments, such as indirect detection or acceler- 
ator searches, are needed in order to break degeneracies between the remaining parameters. 
This illustrates the high degree of complementarity in different searches for SUSY; our results 
show great promise for future direct detection experiments in characterising SUSY together 
with other searches, if it is discovered. 
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